#ifndef _QUADRATURE_F77_DECLARATIONS_H_
#define _QUADRATURE_F77_DECLARATIONS_H_

/*
@brief Raw prototypes of QUADPACK Subroutines
*/

#if defined __cplusplus

extern "C"
{
	/*
	@brief  subroutine DQAG(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
         limit,lenw,last,iwork,work)
		 computation of a definite integral
		 c        standard fortran subroutine
		 c        double precision version
		 c
		 c            f      - double precision
		 c                     function subprogam defining the integrand
		 c                     function f(x). the actual name for f needs to be
		 c                     declared e x t e r n a l in the driver program.
		 c
		 c            a      - double precision
		 c                     lower limit of integration
		 c
		 c            b      - double precision
		 c                     upper limit of integration
		 c
		 c            epsabs - double precision
		 c                     absolute accoracy requested
		 c            epsrel - double precision
		 c                     relative accuracy requested
		 c                     if  epsabs.le.0
		 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
		 c                     the routine will end with ier = 6.
		 c
		 c            key    - integer
		 c                     key for choice of local integration rule
		 c                     a gauss-kronrod pair is used with
		 c                       7 - 15 points if key.lt.2,
		 c                      10 - 21 points if key = 2,
		 c                      15 - 31 points if key = 3,
		 c                      20 - 41 points if key = 4,
		 c                      25 - 51 points if key = 5,
		 c                      30 - 61 points if key.gt.5.
		 c
		 c         on return
		 c            result - double precision
		 c                     approximation to the integral
		 c
		 c            abserr - double precision
		 c                     estimate of the modulus of the absolute error,
		 c                     which should equal or exceed abs(i-result)
		 c
		 c            neval  - integer
		 c                     number of integrand evaluations
		 c
		 c            ier    - integer
		 c                     ier = 0 normal and reliable termination of the
		 c                             routine. it is assumed that the requested
		 c                             accuracy has been achieved.
		 c                     ier.gt.0 abnormal termination of the routine
		 c                             the estimates for result and error are
		 c                             less reliable. it is assumed that the
		 c                             requested accuracy has not been achieved.
		 c                      error messages
		 c                     ier = 1 maximum number of subdivisions allowed
		 c                             has been achieved. one can allow more
		 c                             subdivisions by increasing the value of
		 c                             limit (and taking the according dimension
		 c                             adjustments into account). however, if
		 c                             this yield no improvement it is advised
		 c                             to analyze the integrand in order to
		 c                             determine the integration difficulaties.
		 c                             if the position of a local difficulty can
		 c                             be determined (i.e.singularity,
		 c                             discontinuity within the interval) one
		 c                             will probably gain from splitting up the
		 c                             interval at this point and calling the
		 c                             integrator on the subranges. if possible,
		 c                             an appropriate special-purpose integrator
		 c                             should be used which is designed for
		 c                             handling the type of difficulty involved.
		 c                         = 2 the occurrence of roundoff error is
		 c                             detected, which prevents the requested
		 c                             tolerance from being achieved.
		 c                         = 3 extremely bad integrand behaviour occurs
		 c                             at some points of the integration
		 c                             interval.
		 c                         = 6 the input is invalid, because
		 c                             (epsabs.le.0 and
		 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
		 c                             or limit.lt.1 or lenw.lt.limit*4.
		 c                             result, abserr, neval, last are set
		 c                             to zero.
		 c                             except when lenw is invalid, iwork(1),
		 c                             work(limit*2+1) and work(limit*3+1) are
		 c                             set to zero, work(1) is set to a and
		 c                             work(limit+1) to b.
		 c
		 c         dimensioning parameters
		 c            limit - integer
		 c                    dimensioning parameter for iwork
		 c                    limit determines the maximum number of subintervals
		 c                    in the partition of the given integration interval
		 c                    (a,b), limit.ge.1.
		 c                    if limit.lt.1, the routine will end with ier = 6.
		 c
		 c            lenw  - integer
		 c                    dimensioning parameter for work
		 c                    lenw must be at least limit*4.
		 c                    if lenw.lt.limit*4, the routine will end with
		 c                    ier = 6.
		 c
		 c            last  - integer
		 c                    on return, last equals the number of subintervals
		 c                    produced in the subdiviosion process, which
		 c                    determines the number of significant elements
		 c                    actually in the work arrays.
		 c
		 c         work arrays
		 c            iwork - integer
		 c                    vector of dimension at least limit, the first k
		 c                    elements of which contain pointers to the error
		 c                    estimates over the subintervals, such that
		 c                    work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
		 c                    form a decreasing sequence with k = last if
		 c                    last.le.(limit/2+2), and k = limit+1-last otherwise
		 c
		 c            work  - double precision
		 c                    vector of dimension at least lenw
		 c                    on return
		 c                    work(1), ..., work(last) contain the left end
		 c                    points of the subintervals in the partition of
		 c                     (a,b),
		 c                    work(limit+1), ..., work(limit+last) contain the
		 c                     right end points,
		 c                    work(limit*2+1), ..., work(limit*2+last) contain
		 c                     the integral approximations over the subintervals,
		 c                    work(limit*3+1), ..., work(limit*3+last) contain
		 c                     the error estimates.
	*/
void      DQAG(double(*)(double *), double *, double *, double *, double *, int *, double *,
	double *, int *, int *, int *, int *, int *, int *, double *);

/*
@brief  subroutine DQAGE(f,a,b,epsabs,epsrel,key,limit,result,abserr,
     *   neval,ier,alist,blist,rlist,elist,iord,last)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            key    - integer
	 c                     key for choice of local integration rule
	 c                     a gauss-kronrod pair is used with
	 c                          7 - 15 points if key.lt.2,
	 c                         10 - 21 points if key = 2,
	 c                         15 - 31 points if key = 3,
	 c                         20 - 41 points if key = 4,
	 c                         25 - 51 points if key = 5,
	 c                         30 - 61 points if key.gt.5.
	 c
	 c            limit  - integer
	 c                     gives an upperbound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for result and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value
	 c                             of limit.
	 c                             however, if this yields no improvement it
	 c                             is rather advised to analyze the integrand
	 c                             in order to determine the integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined(e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling the integrator on the
	 c                             subranges. if possible, an appropriate
	 c                             special-purpose integrator should be used
	 c                             which is designed for handling the type of
	 c                             difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             result, abserr, neval, last, rlist(1) ,
	 c                             elist(1) and iord(1) are set to zero.
	 c                             alist(1) and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the left
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            blist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the right
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            rlist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the
	 c                      integral approximations on the subintervals
	 c
	 c            elist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the moduli of the
	 c                      absolute error estimates on the subintervals
	 c
	 c            iord    - integer
	 c                      vector of dimension at least limit, the first k
	 c                      elements of which are pointers to the
	 c                      error estimates over the subintervals,
	 c                      such that elist(iord(1)), ...,
	 c                      elist(iord(k)) form a decreasing sequence,
	 c                      with k = last if last.le.(limit/2+2), and
	 c                      k = limit+1-last otherwise
	 c
	 c            last    - integer
	 c                      number of subintervals actually produced in the
	 c                      subdivision process
*/

void        DQAGE(double(*)(double *), double *, double *, double *, double *, int *, int *, double *,
	double *, int *, int *, double *, double *, double *, double *, int *, int *);

  
/*
"brief   subroutine DQAGI(f,bound,inf,epsabs,epsrel,result,abserr,neval,
     *   ier,limit,lenw,last,iwork,work)
	 integration over infinite intervals
	 c        standard fortran subroutine
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            bound  - double precision
	 c                     finite bound of integration range
	 c                     (has no meaning if interval is doubly-infinite)
	 c
	 c            inf    - integer
	 c                     indicating the kind of integration range involved
	 c                     inf = 1 corresponds to  (bound,+infinity),
	 c                     inf = -1            to  (-infinity,bound),
	 c                     inf = 2             to (-infinity,+infinity).
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for result and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is assumed that the requested tolerance
	 c                             cannot be achieved, and that the returned
	 c                             result is the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                              or limit.lt.1 or leniw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when limit or leniw is
	 c                             invalid, iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first
	 c                    k elements of which contain pointers
	 c                    to the error estimates over the subintervals,
	 c                    such that work(limit*3+iwork(1)),... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2), and
	 c                    k = limit+1-last otherwise
	 c
	 c            work  - double precision
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ...,work(limit*2+last) contain the
	 c                     integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3)
	 c                     contain the error estimates.
*/
void             DQAGI(double(*)(double *), double *, int *, double *, double *, double *, double *,
	int *, int *, int *, int *, int *, int *, double *);

/*
@brief  subroutine DQAGIE(f,bound,inf,epsabs,epsrel,limit,result,abserr,
     *   neval,ier,alist,blist,rlist,elist,iord,last)

	 integration over infinite intervals
	 c standard fortran subroutine
	 c
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            bound  - double precision
	 c                     finite bound of integration range
	 c                     (has no meaning if interval is doubly-infinite)
	 c
	 c            inf    - double precision
	 c                     indicating the kind of integration range involved
	 c                     inf = 1 corresponds to  (bound,+infinity),
	 c                     inf = -1            to  (-infinity,bound),
	 c                     inf = 2             to (-infinity,+infinity).
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for result and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however,if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty can
	 c                             be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is assumed that the requested tolerance
	 c                             cannot be achieved, and that the returned
	 c                             result is the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             result, abserr, neval, last, rlist(1),
	 c                             elist(1) and iord(1) are set to zero.
	 c                             alist(1) and blist(1) are set to 0
	 c                             and 1 respectively.
	 c
	 c            alist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the transformed integration range (0,1).
	 c
	 c            blist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the transformed integration range (0,1).
	 c
	 c            rlist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - double precision
	 c                     vector of dimension at least limit,  the first
	 c                     last elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced
	 c                     in the subdivision process
*/
void           DQAGIE(double(*)(double *), double *, double *, double *, double *, int *, double *,
	double *, int *, int *, double *, double *, double *, double *, int *, int *);

/*
@brief  subroutine DQAGP(f,a,b,npts2,points,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,lenw,last,iwork,work)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            npts2  - integer
	 c                     number equal to two more than the number of
	 c                     user-supplied break points within the integration
	 c                     range, npts.ge.2.
	 c                     if npts2.lt.2, the routine will end with ier = 6.
	 c
	 c            points - double precision
	 c                     vector of dimension npts2, the first (npts2-2)
	 c                     elements of which are the user provided break
	 c                     points. if these points do not constitute an
	 c                     ascending sequence there will be an automatic
	 c                     sorting.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (i.e. singularity,
	 c                             discontinuity within the interval), it
	 c                             should be supplied to the routine as an
	 c                             element of the vector points. if necessary
	 c                             an appropriate special-purpose integrator
	 c                             must be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved, and that
	 c                             the returned result is the best which
	 c                             can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid because
	 c                             npts2.lt.2 or
	 c                             break points are specified outside
	 c                             the integration range or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when leniw or lenw or npts2 is
	 c                             invalid, iwork(1), iwork(limit+1),
	 c                             work(limit*2+1) and work(limit*3+1)
	 c                             are set to zero.
	 c                             work(1) is set to a and work(limit+1)
	 c                             to b (where limit = (leniw-npts2)/2).
	 c
	 c         dimensioning parameters
	 c            leniw - integer
	 c                    dimensioning parameter for iwork
	 c                    leniw determines limit = (leniw-npts2)/2,
	 c                    which is the maximum number of subintervals in the
	 c                    partition of the given integration interval (a,b),
	 c                    leniw.ge.(3*npts2-2).
	 c                    if leniw.lt.(3*npts2-2), the routine will end with
	 c                    ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least leniw*2-npts2.
	 c                    if lenw.lt.leniw*2-npts2, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least leniw. on return,
	 c                    the first k elements of which contain
	 c                    pointers to the error estimates over the
	 c                    subintervals, such that work(limit*3+iwork(1)),...,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2), and
	 c                    k = limit+1-last otherwise
	 c                    iwork(limit+1), ...,iwork(limit+last) contain the
	 c                     subdivision levels of the subintervals, i.e.
	 c                     if (aa,bb) is a subinterval of (p1,p2)
	 c                     where p1 as well as p2 is a user-provided
	 c                     break point or integration limit, then (aa,bb) has
	 c                     level l if abs(bb-aa) = abs(p2-p1)*2**(-l),
	 c                    iwork(limit*2+1), ..., iwork(limit*2+npts2) have
	 c                     no significance for the user,
	 c                    note that limit = (leniw-npts2)/2.
	 c
	 c            work  - double precision
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the corresponding error estimates,
	 c                    work(limit*4+1), ..., work(limit*4+npts2)
	 c                     contain the integration limits and the
	 c                     break points sorted in an ascending sequence.
	 c                    note that limit = (leniw-npts2)/2.
*/
void      DQAGP(double(*)(double *), double *, double *, int *, double *, double *, double *, double *,
	double *, int *, int *, int *, int *, int *, int *, double *);

/*
@brief  subroutine DQAGPE(f,a,b,npts2,points,epsabs,epsrel,limit,result,
     *   abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,
     *   last)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            npts2  - integer
	 c                     number equal to two more than the number of
	 c                     user-supplied break points within the integration
	 c                     range, npts2.ge.2.
	 c                     if npts2.lt.2, the routine will end with ier = 6.
	 c
	 c            points - double precision
	 c                     vector of dimension npts2, the first (npts2-2)
	 c                     elements of which are the user provided break
	 c                     points. if these points do not constitute an
	 c                     ascending sequence there will be an automatic
	 c                     sorting.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.npts2
	 c                     if limit.lt.npts2, the routine will end with
	 c                     ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (i.e. singularity,
	 c                             discontinuity within the interval), it
	 c                             should be supplied to the routine as an
	 c                             element of the vector points. if necessary
	 c                             an appropriate special-purpose integrator
	 c                             must be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be
	 c                             achieved, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid because
	 c                             npts2.lt.2 or
	 c                             break points are specified outside
	 c                             the integration range or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.npts2.
	 c                             result, abserr, neval, last, rlist(1),
	 c                             and elist(1) are set to zero. alist(1) and
	 c                             blist(1) are set to a and b respectively.
	 c
	 c            alist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            blist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            rlist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            pts    - double precision
	 c                     vector of dimension at least npts2, containing the
	 c                     integration limits and the break points of the
	 c                     interval in ascending sequence.
	 c
	 c            level  - integer
	 c                     vector of dimension at least limit, containing the
	 c                     subdivision levels of the subinterval, i.e. if
	 c                     (aa,bb) is a subinterval of (p1,p2) where p1 as
	 c                     well as p2 is a user-provided break point or
	 c                     integration limit, then (aa,bb) has level l if
	 c                     abs(bb-aa) = abs(p2-p1)*2**(-l).
	 c
	 c            ndin   - integer
	 c                     vector of dimension at least npts2, after first
	 c                     integration over the intervals (pts(i)),pts(i+1),
	 c                     i = 0,1, ..., npts2-2, the error estimates over
	 c                     some of the intervals may have been increased
	 c                     artificially, in order to put their subdivision
	 c                     forward. if this happens for the subinterval
	 c                     numbered k, ndin(k) is put to 1, otherwise
	 c                     ndin(k) = 0.
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in the
	 c                     subdivisions process
*/

void        DQAGPE(double(*)(double *), double *, double *, int *, double *, double *, double *,
	int *, double *, double *, int *, int *, double *, double *, double *, double *, double *,
	int *, int *, int *, int *);

/*
@brief   subroutine DQAGS(f,a,b,epsabs,epsrel,result,abserr,neval,ier,
     *   limit,lenw,last,iwork,work)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account. however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be
	 c                             achieved, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28)
	 c                             or limit.lt.1 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero.except when limit or lenw is invalid,
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, detemines the
	 c                    number of significant elements actually in the work
	 c                    arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first k
	 c                    elements of which contain pointers
	 c                    to the error estimates over the subintervals
	 c                    such that work(limit*3+iwork(1)),... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2),
	 c                    and k = limit+1-last otherwise
	 c
	 c            work  - double precision
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end-points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end-points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the error estimates.
*/

void     DQAGS(double(*)(double *), double *, double *, double *, double *, double *, double *,
	int *, int *, int *, int *, int *, int *, double *);

/*
@brief    subroutine DQAGSE(f,a,b,epsabs,epsrel,limit,result,abserr,neval,
     *   ier,alist,blist,rlist,elist,iord,last)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upperbound on the number of subintervals
	 c                     in the partition of (a,b)
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                         = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved, and that the
	 c                             returned result is the best which can be
	 c                             obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             epsabs.le.0 and
	 c                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
	 c                             result, abserr, neval, last, rlist(1),
	 c                             iord(1) and elist(1) are set to zero.
	 c                             alist(1) and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left end points
	 c                     of the subintervals in the partition of the
	 c                     given integration range (a,b)
	 c
	 c            blist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            rlist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in the
	 c                     subdivision process
*/
void             DQAGSE(double(*)(double *), double *, double *, double *, double *, int *, double *,
	double *, int *, int *, double *, double *, double *, double *, int *, int *);

/*
@brief   subroutine DQAWC(f,a,b,c,epsabs,epsrel,result,abserr,neval,ier,
     *   limit,lenw,last,iwork,work)

	 computation of a cauchy principal value
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     under limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            c      - parameter in the weight function, c.ne.a, c.ne.b.
	 c                     if c = a or c = b, the routine will end with
	 c                     ier = 6 .
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate or the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty
	 c                             can be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling
	 c                             appropriate integrators on the subranges.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             c = a or c = b or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.1 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when lenw or limit is invalid,
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c           lenw   - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end with
	 c                    ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first k
	 c                    elements of which contain pointers
	 c                    to the error estimates over the subintervals,
	 c                    such that work(limit*3+iwork(1)), ... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2),
	 c                    and k = limit+1-last otherwise
	 c
	 c            work  - double precision
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the error estimates.
*/

void              DQAWC(double(*)(double *), double *, double *, double *, double *, double *, double *,
	double *, int *, int *, int *, int *, int *, int *, double *);

/*
@brief     subroutine DQAWCE(f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,
     *   ier,alist,blist,rlist,elist,iord,last)

	 computation of a cauchy principal value
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            c      - double precision
	 c                     parameter in the weight function, c.ne.a, c.ne.b
	 c                     if c = a or c = b, the routine will end with
	 c                     ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of
	 c                             limit. however, if this yields no
	 c                             improvement it is advised to analyze the
	 c                             the integrand, in order to determine the
	 c                             the integration difficulties. if the
	 c                             position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling
	 c                             appropriate integrators on the subranges.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some interior points of
	 c                             the integration interval.
	 c                         = 6 the input is invalid, because
	 c                             c = a or c = b or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.1.
	 c                             result, abserr, neval, rlist(1), elist(1),
	 c                             iord(1) and last are set to zero. alist(1)
	 c                             and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the left
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            blist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the right
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            rlist   - double precision
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the integral
	 c                      approximations on the subintervals
	 c
	 c            elist   - double precision
	 c                      vector of dimension limit, the first  last
	 c                      elements of which are the moduli of the absolute
	 c                      error estimates on the subintervals
	 c
	 c            iord    - integer
	 c                      vector of dimension at least limit, the first k
	 c                      elements of which are pointers to the error
	 c                      estimates over the subintervals, so that
	 c                      elist(iord(1)), ..., elist(iord(k)) with k = last
	 c                      if last.le.(limit/2+2), and k = limit+1-last
	 c                      otherwise, form a decreasing sequence
	 c
	 c            last    - integer
	 c                      number of subintervals actually produced in
	 c                      the subdivision process
*/

void               DQAWCE(double(*)(double *), double *, double *, double *, double *, double *, int *,
	double *, double *, int *, int *, double *, double *, double *, double *, int *, int *);

/*
@brief        subroutine DQAWF(f,a,omega,integr,epsabs,result,abserr,neval,ier,
     *   limlst,lst,leniw,maxp1,lenw,iwork,work)

	 computation of fourier integrals
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            omega  - double precision
	 c                     parameter in the integrand weight function
	 c
	 c            integr - integer
	 c                     indicates which of the weight functions is used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1.and.integr.ne.2, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested, epsabs.gt.0.
	 c                     if epsabs.le.0, the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                    if omega.ne.0
	 c                     ier = 1 maximum number of cycles allowed
	 c                             has been achieved, i.e. of subintervals
	 c                             (a+(k-1)c,a+kc) where
	 c                             c = (2*int(abs(omega))+1)*pi/abs(omega),
	 c                             for k = 1, 2, ..., lst.
	 c                             one can allow more cycles by increasing
	 c                             the value of limlst (and taking the
	 c                             according dimension adjustments into
	 c                             account). examine the array iwork which
	 c                             contains the error flags on the cycles, in
	 c                             order to look for eventual local
	 c                             integration difficulties.
	 c                             if the position of a local difficulty
	 c                             can be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling
	 c                             appropriate integrators on the subranges.
	 c                         = 4 the extrapolation table constructed for
	 c                             convergence accelaration of the series
	 c                             formed by the integral contributions over
	 c                             the cycles, does not converge to within
	 c                             the requested accuracy.
	 c                             as in the case of ier = 1, it is advised
	 c                             to examine the array iwork which contains
	 c                             the error flags on the cycles.
	 c                         = 6 the input is invalid because
	 c                             (integr.ne.1 and integr.ne.2) or
	 c                              epsabs.le.0 or limlst.lt.1 or
	 c                              leniw.lt.(limlst+2) or maxp1.lt.1 or
	 c                              lenw.lt.(leniw*2+maxp1*25).
	 c                              result, abserr, neval, lst are set to
	 c                              zero.
	 c                         = 7 bad integrand behaviour occurs within
	 c                             one or more of the cycles. location and
	 c                             type of the difficulty involved can be
	 c                             determined from the first lst elements of
	 c                             vector iwork.  here lst is the number of
	 c                             cycles actually needed (see below).
	 c                             iwork(k) = 1 the maximum number of
	 c                                          subdivisions (=(leniw-limlst)
	 c                                          /2) has been achieved on the
	 c                                          k th cycle.
	 c                                      = 2 occurrence of roundoff error
	 c                                          is detected and prevents the
	 c                                          tolerance imposed on the k th
	 c                                          cycle, from being achieved
	 c                                          on this cycle.
	 c                                      = 3 extremely bad integrand
	 c                                          behaviour occurs at some
	 c                                          points of the k th cycle.
	 c                                      = 4 the integration procedure
	 c                                          over the k th cycle does
	 c                                          not converge (to within the
	 c                                          required accuracy) due to
	 c                                          roundoff in the extrapolation
	 c                                          procedure invoked on this
	 c                                          cycle. it is assumed that the
	 c                                          result on this interval is
	 c                                          the best which can be
	 c                                          obtained.
	 c                                      = 5 the integral over the k th
	 c                                          cycle is probably divergent
	 c                                          or slowly convergent. it must
	 c                                          be noted that divergence can
	 c                                          occur with any other value of
	 c                                          iwork(k).
	 c                    if omega = 0 and integr = 1,
	 c                    the integral is calculated by means of dqagie,
	 c                    and ier = iwork(1) (with meaning as described
	 c                    for iwork(k),k = 1).
	 c
	 c         dimensioning parameters
	 c            limlst - integer
	 c                     limlst gives an upper bound on the number of
	 c                     cycles, limlst.ge.3.
	 c                     if limlst.lt.3, the routine will end with ier = 6.
	 c
	 c            lst    - integer
	 c                     on return, lst indicates the number of cycles
	 c                     actually needed for the integration.
	 c                     if omega = 0, then lst is set to 1.
	 c
	 c            leniw  - integer
	 c                     dimensioning parameter for iwork. on entry,
	 c                     (leniw-limlst)/2 equals the maximum number of
	 c                     subintervals allowed in the partition of each
	 c                     cycle, leniw.ge.(limlst+2).
	 c                     if leniw.lt.(limlst+2), the routine will end with
	 c                     ier = 6.
	 c
	 c            maxp1  - integer
	 c                     maxp1 gives an upper bound on the number of
	 c                     chebyshev moments which can be stored, i.e. for
	 c                     the intervals of lengths abs(b-a)*2**(-l),
	 c                     l = 0,1, ..., maxp1-2, maxp1.ge.1.
	 c                     if maxp1.lt.1, the routine will end with ier = 6.
	 c            lenw   - integer
	 c                     dimensioning parameter for work
	 c                     lenw must be at least leniw*2+maxp1*25.
	 c                     if lenw.lt.(leniw*2+maxp1*25), the routine will
	 c                     end with ier = 6.
	 c
	 c         work arrays
	 c            iwork  - integer
	 c                     vector of dimension at least leniw
	 c                     on return, iwork(k) for k = 1, 2, ..., lst
	 c                     contain the error flags on the cycles.
	 c
	 c            work   - double precision
	 c                     vector of dimension at least
	 c                     on return,
	 c                     work(1), ..., work(lst) contain the integral
	 c                      approximations over the cycles,
	 c                     work(limlst+1), ..., work(limlst+lst) contain
	 c                      the error extimates over the cycles.
	 c                     further elements of work have no specific
	 c                     meaning for the user.
*/

void          DQAWF(double(*)(double *), double *, double *, int *, double *, double *, double *,
	int *, int *, int *, int *, int *, int *, int *, int *, double *);

/*
brief     subroutine DQAWFE(f,a,omega,integr,epsabs,limlst,limit,maxp1,
     *   result,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,
     *   rlist,elist,iord,nnlog,chebmo)

	 computation of fourier integrals
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to
	 c                     be declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            omega  - double precision
	 c                     parameter in the weight function
	 c
	 c            integr - integer
	 c                     indicates which weight function is used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1.and.integr.ne.2, the routine will
	 c                     end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested, epsabs.gt.0
	 c                     if epsabs.le.0, the routine will end with ier = 6.
	 c
	 c            limlst - integer
	 c                     limlst gives an upper bound on the number of
	 c                     cycles, limlst.ge.1.
	 c                     if limlst.lt.3, the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     allowed in the partition of each cycle, limit.ge.1
	 c                     each cycle, limit.ge.1.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of
	 c                     chebyshev moments which can be stored, i.e.
	 c                     for the intervals of lengths abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral x
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - ier = 0 normal and reliable termination of
	 c                             the routine. it is assumed that the
	 c                             requested accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for integral and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                    if omega.ne.0
	 c                     ier = 1 maximum number of  cycles  allowed
	 c                             has been achieved., i.e. of subintervals
	 c                             (a+(k-1)c,a+kc) where
	 c                             c = (2*int(abs(omega))+1)*pi/abs(omega),
	 c                             for k = 1, 2, ..., lst.
	 c                             one can allow more cycles by increasing
	 c                             the value of limlst (and taking the
	 c                             according dimension adjustments into
	 c                             account).
	 c                             examine the array iwork which contains
	 c                             the error flags on the cycles, in order to
	 c                             look for eventual local integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined (e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling appropriate integrators on
	 c                             the subranges.
	 c                         = 4 the extrapolation table constructed for
	 c                             convergence acceleration of the series
	 c                             formed by the integral contributions over
	 c                             the cycles, does not converge to within
	 c                             the requested accuracy. as in the case of
	 c                             ier = 1, it is advised to examine the
	 c                             array iwork which contains the error
	 c                             flags on the cycles.
	 c                         = 6 the input is invalid because
	 c                             (integr.ne.1 and integr.ne.2) or
	 c                              epsabs.le.0 or limlst.lt.3.
	 c                              result, abserr, neval, lst are set
	 c                              to zero.
	 c                         = 7 bad integrand behaviour occurs within one
	 c                             or more of the cycles. location and type
	 c                             of the difficulty involved can be
	 c                             determined from the vector ierlst. here
	 c                             lst is the number of cycles actually
	 c                             needed (see below).
	 c                             ierlst(k) = 1 the maximum number of
	 c                                           subdivisions (= limit) has
	 c                                           been achieved on the k th
	 c                                           cycle.
	 c                                       = 2 occurrence of roundoff error
	 c                                           is detected and prevents the
	 c                                           tolerance imposed on the
	 c                                           k th cycle, from being
	 c                                           achieved.
	 c                                       = 3 extremely bad integrand
	 c                                           behaviour occurs at some
	 c                                           points of the k th cycle.
	 c                                       = 4 the integration procedure
	 c                                           over the k th cycle does
	 c                                           not converge (to within the
	 c                                           required accuracy) due to
	 c                                           roundoff in the
	 c                                           extrapolation procedure
	 c                                           invoked on this cycle. it
	 c                                           is assumed that the result
	 c                                           on this interval is the
	 c                                           best which can be obtained.
	 c                                       = 5 the integral over the k th
	 c                                           cycle is probably divergent
	 c                                           or slowly convergent. it
	 c                                           must be noted that
	 c                                           divergence can occur with
	 c                                           any other value of
	 c                                           ierlst(k).
	 c                    if omega = 0 and integr = 1,
	 c                    the integral is calculated by means of dqagie
	 c                    and ier = ierlst(1) (with meaning as described
	 c                    for ierlst(k), k = 1).
	 c
	 c            rslst  - double precision
	 c                     vector of dimension at least limlst
	 c                     rslst(k) contains the integral contribution
	 c                     over the interval (a+(k-1)c,a+kc) where
	 c                     c = (2*int(abs(omega))+1)*pi/abs(omega),
	 c                     k = 1, 2, ..., lst.
	 c                     note that, if omega = 0, rslst(1) contains
	 c                     the value of the integral over (a,infinity).
	 c
	 c            erlst  - double precision
	 c                     vector of dimension at least limlst
	 c                     erlst(k) contains the error estimate corresponding
	 c                     with rslst(k).
	 c
	 c            ierlst - integer
	 c                     vector of dimension at least limlst
	 c                     ierlst(k) contains the error flag corresponding
	 c                     with rslst(k). for the meaning of the local error
	 c                     flags see description of output parameter ier.
	 c
	 c            lst    - integer
	 c                     number of subintervals needed for the integration
	 c                     if omega = 0 then lst is set to 1.
	 c
	 c            alist, blist, rlist, elist - double precision
	 c                     vector of dimension at least limit,
	 c
	 c            iord, nnlog - integer
	 c                     vector of dimension at least limit, providing
	 c                     space for the quantities needed in the subdivision
	 c                     process of each cycle
	 c
	 c            chebmo - double precision
	 c                     array of dimension at least (maxp1,25), providing
	 c                     space for the chebyshev moments needed within the
	 c                     cycles
*/

void            DQAWFE(double(*)(double *), double *, double *, int *, double *, int *, int *, int *,
	double *, double *, int *, int *, double *, double *, int *, int *, double *, double *,
	double *, double *, int *, int *, double *);

/*
@brief      subroutine DQAWO(f,a,b,omega,integr,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,maxp1,lenw,last,iwork,work)

	 computation of oscillatory integrals
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the function
	 c                     f(x).  the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            omega  - double precision
	 c                     parameter in the integrand weight function
	 c
	 c            integr - integer
	 c                     indicates which of the weight functions is used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1.and.integr.ne.2, the routine will
	 c                     end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if epsabs.le.0 and
	 c                     epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of  integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             (= leniw/2) has been achieved. one can
	 c                             allow more subdivisions by increasing the
	 c                             value of leniw (and taking the according
	 c                             dimension adjustments into account).
	 c                             however, if this yields no improvement it
	 c                             is advised to analyze the integrand in
	 c                             order to determine the integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined (e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling the integrator on the
	 c                             subranges. if possible, an appropriate
	 c                             special-purpose integrator should be used
	 c                             which is designed for handling the type of
	 c                             difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some interior points of the
	 c                             integration interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be achieved
	 c                             due to roundoff in the extrapolation
	 c                             table, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or (integr.ne.1 and integr.ne.2),
	 c                             or leniw.lt.2 or maxp1.lt.1 or
	 c                             lenw.lt.leniw*2+maxp1*25.
	 c                             result, abserr, neval, last are set to
	 c                             zero. except when leniw, maxp1 or lenw are
	 c                             invalid, work(limit*2+1), work(limit*3+1),
	 c                             iwork(1), iwork(limit+1) are set to zero,
	 c                             work(1) is set to a and work(limit+1) to
	 c                             b.
	 c
	 c         dimensioning parameters
	 c            leniw  - integer
	 c                     dimensioning parameter for iwork.
	 c                     leniw/2 equals the maximum number of subintervals
	 c                     allowed in the partition of the given integration
	 c                     interval (a,b), leniw.ge.2.
	 c                     if leniw.lt.2, the routine will end with ier = 6.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of chebyshev
	 c                     moments which can be stored, i.e. for the
	 c                     intervals of lengths abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1
	 c                     if maxp1.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw   - integer
	 c                     dimensioning parameter for work
	 c                     lenw must be at least leniw*2+maxp1*25.
	 c                     if lenw.lt.(leniw*2+maxp1*25), the routine will
	 c                     end with ier = 6.
	 c
	 c            last   - integer
	 c                     on return, last equals the number of subintervals
	 c                     produced in the subdivision process, which
	 c                     determines the number of significant elements
	 c                     actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork  - integer
	 c                     vector of dimension at least leniw
	 c                     on return, the first k elements of which contain
	 c                     pointers to the error estimates over the
	 c                     subintervals, such that work(limit*3+iwork(1)), ..
	 c                     work(limit*3+iwork(k)) form a decreasing
	 c                     sequence, with limit = lenw/2 , and k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise.
	 c                     furthermore, iwork(limit+1), ..., iwork(limit+
	 c                     last) indicate the subdivision levels of the
	 c                     subintervals, such that iwork(limit+i) = l means
	 c                     that the subinterval numbered i is of length
	 c                     abs(b-a)*2**(1-l).
	 c
	 c            work   - double precision
	 c                     vector of dimension at least lenw
	 c                     on return
	 c                     work(1), ..., work(last) contain the left
	 c                      end points of the subintervals in the
	 c                      partition of (a,b),
	 c                     work(limit+1), ..., work(limit+last) contain
	 c                      the right end points,
	 c                     work(limit*2+1), ..., work(limit*2+last) contain
	 c                      the integral approximations over the
	 c                      subintervals,
	 c                     work(limit*3+1), ..., work(limit*3+last)
	 c                      contain the error estimates.
	 c                     work(limit*4+1), ..., work(limit*4+maxp1*25)
	 c                      provide space for storing the chebyshev moments.
	 c                     note that limit = lenw/2.
	 c
*/

void             DQAWO(double(*)(double *), double *, double *, double *, int *, double *, double *, double *, double *,
	int *, int *, int *, int *, int *, int *, int *, double *);

/*
@brief       subroutine  DQAWOE(f,a,b,omega,integr,epsabs,epsrel,limit,icall,
     *  maxp1,result,abserr,neval,ier,last,alist,blist,rlist,elist,iord,
     *   nnlog,momcom,chebmo)

	 computation of oscillatory integrals
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration
	 c
	 c            omega  - double precision
	 c                     parameter in the integrand weight function
	 c
	 c            integr - integer
	 c                     indicates which of the weight functions is to be
	 c                     used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1 and integr.ne.2, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subdivisions
	 c                     in the partition of (a,b), limit.ge.1.
	 c
	 c            icall  - integer
	 c                     if dqawoe is to be used only once, icall must
	 c                     be set to 1.  assume that during this call, the
	 c                     chebyshev moments (for clenshaw-curtis integration
	 c                     of degree 24) have been computed for intervals of
	 c                     lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1.
	 c                     if icall.gt.1 this means that dqawoe has been
	 c                     called twice or more on intervals of the same
	 c                     length abs(b-a). the chebyshev moments already
	 c                     computed are then re-used in subsequent calls.
	 c                     if icall.lt.1, the routine will end with ier = 6.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of chebyshev
	 c                     moments which can be stored, i.e. for the
	 c                     intervals of lenghts abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1.
	 c                     if maxp1.lt.1, the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the
	 c                             requested accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand, in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty can
	 c                             be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved due to
	 c                             roundoff in the extrapolation table,
	 c                             and that the returned result is the
	 c                             best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or (integr.ne.1 and integr.ne.2) or
	 c                             icall.lt.1 or maxp1.lt.1.
	 c                             result, abserr, neval, last, rlist(1),
	 c                             elist(1), iord(1) and nnlog(1) are set
	 c                             to zero. alist(1) and blist(1) are set
	 c                             to a and b respectively.
	 c
	 c            last  -  integer
	 c                     on return, last equals the number of
	 c                     subintervals produces in the subdivision
	 c                     process, which determines the number of
	 c                     significant elements actually in the
	 c                     work arrays.
	 c            alist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            blist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            rlist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the error
	 c                     estimates over the subintervals,
	 c                     such that elist(iord(1)), ...,
	 c                     elist(iord(k)) form a decreasing sequence, with
	 c                     k = last if last.le.(limit/2+2), and
	 c                     k = limit+1-last otherwise.
	 c
	 c            nnlog  - integer
	 c                     vector of dimension at least limit, containing the
	 c                     subdivision levels of the subintervals, i.e.
	 c                     iwork(i) = l means that the subinterval
	 c                     numbered i is of length abs(b-a)*2**(1-l)
	 c
	 c         on entry and return
	 c            momcom - integer
	 c                     indicating that the chebyshev moments
	 c                     have been computed for intervals of lengths
	 c                     (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1,
	 c                     momcom.lt.maxp1
	 c
	 c            chebmo - double precision
	 c                     array of dimension (maxp1,25) containing the
	 c                     chebyshev moments
*/

void              DQAWOE(double(*)(double *), double *, double *, double *, int *, double *, double *, int *,
	int *, int *, double *, double *, int *, int *, int *, double *, double *, double *,
	double *, int *, int *, int *, double *);

/*
@brief       subroutine DQAWS(f,a,b,alfa,beta,integr,epsabs,epsrel,result,
     *   abserr,neval,ier,limit,lenw,last,iwork,work)

	 integration of functions having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration, b.gt.a
	 c                     if b.le.a, the routine will end with ier = 6.
	 c
	 c            alfa   - double precision
	 c                     parameter in the integrand function, alfa.gt.(-1)
	 c                     if alfa.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            beta   - double precision
	 c                     parameter in the integrand function, beta.gt.(-1)
	 c                     if beta.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            integr - integer
	 c                     indicates which weight function is to be used
	 c                     = 1  (x-a)**alfa*(b-x)**beta
	 c                     = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                     = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                     = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
	 c                     if integr.lt.1 or integr.gt.4, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal tep
	 c                       ination of the routine
	 c                             the estimates for the integral and error
	 c                             are less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand, in order to
	 c                             determine the integration difficulties
	 c                             which prevent the requested tolerance from
	 c                             being achieved. in case of a jump
	 c                             discontinuity or a local singularity
	 c                             of algebraico-logarithmic type at one or
	 c                             more interior points of the integration
	 c                             range, one should proceed by splitting up
	 c                             the interval at these points and calling
	 c                             the integrator on the subranges.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             b.le.a or alfa.le.(-1) or beta.le.(-1) or
	 c                             or integr.lt.1 or integr.gt.4 or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.2 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. except when lenw or limit is invalid
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit  - integer
	 c                     dimensioning parameter for iwork
	 c                     limit determines the maximum number of
	 c                     subintervals in the partition of the given
	 c                     integration interval (a,b), limit.ge.2.
	 c                     if limit.lt.2, the routine will end with ier = 6.
	 c
	 c            lenw   - integer
	 c                     dimensioning parameter for work
	 c                     lenw must be at least limit*4.
	 c                     if lenw.lt.limit*4, the routine will end
	 c                     with ier = 6.
	 c
	 c            last   - integer
	 c                     on return, last equals the number of
	 c                     subintervals produced in the subdivision process,
	 c                     which determines the significant number of
	 c                     elements actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork  - integer
	 c                     vector of dimension limit, the first k
	 c                     elements of which contain pointers
	 c                     to the error estimates over the subintervals,
	 c                     such that work(limit*3+iwork(1)), ...,
	 c                     work(limit*3+iwork(k)) form a decreasing
	 c                     sequence with k = last if last.le.(limit/2+2),
	 c                     and k = limit+1-last otherwise
	 c
	 c            work   - double precision
	 c                     vector of dimension lenw
	 c                     on return
	 c                     work(1), ..., work(last) contain the left
	 c                      end points of the subintervals in the
	 c                      partition of (a,b),
	 c                     work(limit+1), ..., work(limit+last) contain
	 c                      the right end points,
	 c                     work(limit*2+1), ..., work(limit*2+last)
	 c                      contain the integral approximations over
	 c                      the subintervals,
	 c                     work(limit*3+1), ..., work(limit*3+last)
	 c                      contain the error estimates.
*/

void              DQAWS(double(*)(double *), double *, double *, double *, double *, int *, double *,
	double *, double *, double *, int *, int *, int *, int *, int *, int *, double *);

/*
brief  subroutine DQAWSE(f,a,b,alfa,beta,integr,epsabs,epsrel,limit,
     *   result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)

	 integration of functions having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c            f      - double precision
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - double precision
	 c                     lower limit of integration
	 c
	 c            b      - double precision
	 c                     upper limit of integration, b.gt.a
	 c                     if b.le.a, the routine will end with ier = 6.
	 c
	 c            alfa   - double precision
	 c                     parameter in the weight function, alfa.gt.(-1)
	 c                     if alfa.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            beta   - double precision
	 c                     parameter in the weight function, beta.gt.(-1)
	 c                     if beta.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            integr - integer
	 c                     indicates which weight function is to be used
	 c                     = 1  (x-a)**alfa*(b-x)**beta
	 c                     = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                     = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                     = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
	 c                     if integr.lt.1 or integr.gt.4, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - double precision
	 c                     absolute accuracy requested
	 c            epsrel - double precision
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.2
	 c                     if limit.lt.2, the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - double precision
	 c                     approximation to the integral
	 c
	 c            abserr - double precision
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for the integral and error
	 c                             are less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                         = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit. however, if this yields no
	 c                             improvement, it is advised to analyze the
	 c                             integrand in order to determine the
	 c                             integration difficulties which prevent the
	 c                             requested tolerance from being achieved.
	 c                             in case of a jump discontinuity or a local
	 c                             singularity of algebraico-logarithmic type
	 c                             at one or more interior points of the
	 c                             integration range, one should proceed by
	 c                             splitting up the interval at these
	 c                             points and calling the integrator on the
	 c                             subranges.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             b.le.a or alfa.le.(-1) or beta.le.(-1), or
	 c                             integr.lt.1 or integr.gt.4, or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             or limit.lt.2.
	 c                             result, abserr, neval, rlist(1), elist(1),
	 c                             iord(1) and last are set to zero. alist(1)
	 c                             and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            blist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            rlist  - double precision
	 c                     vector of dimension at least limit,the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - double precision
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     of which are pointers to the error
	 c                     estimates over the subintervals, so that
	 c                     elist(iord(1)), ..., elist(iord(k)) with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise form a decreasing sequence
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in
	 c                     the subdivision process
*/

void               DQAWSE(double(*)(double *), double *, double *, double *, double *, int *, double *,
	double *, int *, double *, double *, int *, int *, double *, double *, double *, double *, int *, int *);

/*
@brief    subroutine DQC25C(f,a,b,c,result,abserr,krul,neval)

    integration rules for the computation of cauchy
c        principal value integrals
c        standard fortran subroutine
c        double precision version
c
c        parameters
c           f      - double precision
c                    function subprogram defining the integrand function
c                    f(x). the actual name for f needs to be declared
c                    e x t e r n a l  in the driver program.
c
c           a      - double precision
c                    left end point of the integration interval
c
c           b      - double precision
c                    right end point of the integration interval, b.gt.a
c
c           c      - double precision
c                    parameter in the weight function
c
c           result - double precision
c                    approximation to the integral
c                    result is computed by using a generalized
c                    clenshaw-curtis method if c lies within ten percent
c                    of the integration interval. in the other case the
c                    15-point kronrod rule obtained by optimal addition
c                    of abscissae to the 7-point gauss rule, is applied.
c
c           abserr - double precision
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed abs(i-result)
c
c           krul   - integer
c                    key which is decreased by 1 if the 15-point
c                    gauss-kronrod scheme has been used
c
c           neval  - integer
c                    number of integrand evaluations
c
*/

void        DQC25C(double(*)(double *), double *, double *, double *, double *, double *, int *, int *);

/*
@brief      subroutine DQC25F(f,a,b,omega,integr,nrmom,maxp1,ksave,result,
     *   abserr,neval,resabs,resasc,momcom,chebmo)

	 integration rules for functions with cos or sin factor
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c         on entry
	 c           f      - double precision
	 c                    function subprogram defining the integrand
	 c                    function f(x). the actual name for f needs to
	 c                    be declared e x t e r n a l in the calling program.
	 c
	 c           a      - double precision
	 c                    lower limit of integration
	 c
	 c           b      - double precision
	 c                    upper limit of integration
	 c
	 c           omega  - double precision
	 c                    parameter in the weight function
	 c
	 c           integr - integer
	 c                    indicates which weight function is to be used
	 c                       integr = 1   w(x) = cos(omega*x)
	 c                       integr = 2   w(x) = sin(omega*x)
	 c
	 c           nrmom  - integer
	 c                    the length of interval (a,b) is equal to the length
	 c                    of the original integration interval divided by
	 c                    2**nrmom (we suppose that the routine is used in an
	 c                    adaptive integration process, otherwise set
	 c                    nrmom = 0). nrmom must be zero at the first call.
	 c
	 c           maxp1  - integer
	 c                    gives an upper bound on the number of chebyshev
	 c                    moments which can be stored, i.e. for the
	 c                    intervals of lengths abs(bb-aa)*2**(-l),
	 c                    l = 0,1,2, ..., maxp1-2.
	 c
	 c           ksave  - integer
	 c                    key which is one when the moments for the
	 c                    current interval have been computed
	 c
	 c         on return
	 c           result - double precision
	 c                    approximation to the integral i
	 c
	 c           abserr - double precision
	 c                    estimate of the modulus of the absolute
	 c                    error, which should equal or exceed abs(i-result)
	 c
	 c           neval  - integer
	 c                    number of integrand evaluations
	 c
	 c           resabs - double precision
	 c                    approximation to the integral j
	 c
	 c           resasc - double precision
	 c                    approximation to the integral of abs(f-i/(b-a))
	 c
	 c         on entry and return
	 c           momcom - integer
	 c                    for each interval length we need to compute the
	 c                    chebyshev moments. momcom counts the number of
	 c                    intervals for which these moments have already been
	 c                    computed. if nrmom.lt.momcom or ksave = 1, the
	 c                    chebyshev moments for the interval (a,b) have
	 c                    already been computed and stored, otherwise we
	 c                    compute them and we increase momcom.
	 c
	 c           chebmo - double precision
	 c                    array of dimension at least (maxp1,25) containing
	 c                    the modified chebyshev moments for the first momcom
	 c                    momcom interval lengths
*/

void             DQC25F(double(*)(double *), double *, double *, double *, int *, int *, int *,
	int *, double *, double *, int *, double *, double *, int *, double *);

/*
@brief     subroutine DQC25S(f,a,b,bl,br,alfa,beta,ri,rj,rg,rh,result,
     *   abserr,resasc,integr,nev)

	 integration rules for integrands having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        double precision version
	 c
	 c        parameters
	 c           f      - double precision
	 c                    function subprogram defining the integrand
	 c                    f(x). the actual name for f needs to be declared
	 c                    e x t e r n a l  in the driver program.
	 c
	 c           a      - double precision
	 c                    left end point of the original interval
	 c
	 c           b      - double precision
	 c                    right end point of the original interval, b.gt.a
	 c
	 c           bl     - double precision
	 c                    lower limit of integration, bl.ge.a
	 c
	 c           br     - double precision
	 c                    upper limit of integration, br.le.b
	 c
	 c           alfa   - double precision
	 c                    parameter in the weight function
	 c
	 c           beta   - double precision
	 c                    parameter in the weight function
	 c
	 c           ri,rj,rg,rh - double precision
	 c                    modified chebyshev moments for the application
	 c                    of the generalized clenshaw-curtis
	 c                    method (computed in subroutine dqmomo)
	 c
	 c           result - double precision
	 c                    approximation to the integral
	 c                    result is computed by using a generalized
	 c                    clenshaw-curtis method if b1 = a or br = b.
	 c                    in all other cases the 15-point kronrod
	 c                    rule is applied, obtained by optimal addition of
	 c                    abscissae to the 7-point gauss rule.
	 c
	 c           abserr - double precision
	 c                    estimate of the modulus of the absolute error,
	 c                    which should equal or exceed abs(i-result)
	 c
	 c           resasc - double precision
	 c                    approximation to the integral of abs(f*w-i/(b-a))
	 c
	 c           integr - integer
	 c                    which determines the weight function
	 c                    = 1   w(x) = (x-a)**alfa*(b-x)**beta
	 c                    = 2   w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                    = 3   w(x) = (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                    = 4   w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)*
	 c                                 log(b-x)
	 c
	 c           nev    - integer
	 c                    number of integrand evaluations
*/

void            DQC25S(double(*)(double *), double *, double *, double *, double *, double *, double *,
	double *, double *, double *, double *, double *, double *, double *, int *, int *);

/*
@brief     subroutine DQCHEB(x,fval,cheb12,cheb24)

    chebyshev series expansion
c        standard fortran subroutine
c        double precision version
c
c        parameters
c          on entry
c           x      - double precision
c                    vector of dimension 11 containing the
c                    values cos(k*pi/24), k = 1, ..., 11
c
c           fval   - double precision
c                    vector of dimension 25 containing the
c                    function values at the points
c                    (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24,
c                    where (a,b) is the approximation interval.
c                    fval(1) and fval(25) are divided by two
c                    (these values are destroyed at output).
c
c          on return
c           cheb12 - double precision
c                    vector of dimension 13 containing the
c                    chebyshev coefficients for degree 12
c
c           cheb24 - double precision
c                    vector of dimension 25 containing the
c                    chebyshev coefficients for degree 24
*/

void        DQCHEB(double *, double *, double *, double *);

/*
@brief   subroutine DQELG(n,epstab,result,abserr,res3la,nres)

epsilon algorithm
c           standard fortran subroutine
c           double precision version
c
c           parameters
c              n      - integer
c                       epstab(n) contains the new element in the
c                       first column of the epsilon table.
c
c              epstab - double precision
c                       vector of dimension 52 containing the elements
c                       of the two lower diagonals of the triangular
c                       epsilon table. the elements are numbered
c                       starting at the right-hand corner of the
c                       triangle.
c
c              result - double precision
c                       resulting approximation to the integral
c
c              abserr - double precision
c                       estimate of the absolute error computed from
c                       result and the 3 previous results
c
c              res3la - double precision
c                       vector of dimension 3 containing the last 3
c                       results
c
c              nres   - integer
c                       number of calls to the routine
c                       (should be zero at first call)
*/

void          DQELG(int *, double *, double *, double *, double *, int *);

/*
@brief   subroutine DQK15(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the7-point gauss rule(resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void            DQK15(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief    subroutine DQK15I(f,boun,inf,a,b,result,abserr,resabs,resasc)

integration rule
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       fuction subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              boun   - double precision
c                       finite bound of original integration
c                       range (set to zero if inf = +2)
c
c              inf    - integer
c                       if inf = -1, the original interval is
c                                   (-infinity,bound),
c                       if inf = +1, the original interval is
c                                   (bound,+infinity),
c                       if inf = +2, the original interval is
c                                   (-infinity,+infinity) and
c                       the integral is computed as the sum of two
c                       integrals, one over (-infinity,0) and one over
c                       (0,+infinity).
c
c              a      - double precision
c                       lower limit for integration over subrange
c                       of (0,1)
c
c              b      - double precision
c                       upper limit for integration over subrange
c                       of (0,1)
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule(resk) obtained by optimal addition
c                       of abscissae to the 7-point gauss rule(resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should equal or exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of
c                       abs((transformed integrand)-i/(b-a)) over (a,b)
c
*/

void           DQK15I(double(*)(double *), double *, int *, double *, double *, double *, double *, double *, double *);

/*
@brief    subroutine DQK15W(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,
     *   resabs,resasc)

	 integration rules
	 c           standard fortran subroutine
	 c           double precision version
	 c
	 c           parameters
	 c             on entry
	 c              f      - double precision
	 c                       function subprogram defining the integrand
	 c                       function f(x). the actual name for f needs to be
	 c                       declared e x t e r n a l in the driver program.
	 c
	 c              w      - double precision
	 c                       function subprogram defining the integrand
	 c                       weight function w(x). the actual name for w
	 c                       needs to be declared e x t e r n a l in the
	 c                       calling program.
	 c
	 c              p1, p2, p3, p4 - double precision
	 c                       parameters in the weight function
	 c
	 c              kp     - integer
	 c                       key for indicating the type of weight function
	 c
	 c              a      - double precision
	 c                       lower limit of integration
	 c
	 c              b      - double precision
	 c                       upper limit of integration
	 c
	 c            on return
	 c              result - double precision
	 c                       approximation to the integral i
	 c                       result is computed by applying the 15-point
	 c                       kronrod rule (resk) obtained by optimal addition
	 c                       of abscissae to the 7-point gauss rule (resg).
	 c
	 c              abserr - double precision
	 c                       estimate of the modulus of the absolute error,
	 c                       which should equal or exceed abs(i-result)
	 c
	 c              resabs - double precision
	 c                       approximation to the integral of abs(f)
	 c
	 c              resasc - double precision
	 c                       approximation to the integral of abs(f-i/(b-a))
	 c
*/

void              DQK15W(double(*)(double *), double *, double *, double *, double *, double *,
	int *, double *, double *, double *, double *, double *, double *);

/*
@brief   subroutine DQK21(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the driver program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 21-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the 10-point gauss rule (resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void           DQK21(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief    subroutine DQK31(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 31-point
c                       gauss-kronrod rule (resk), obtained by optimal
c                       addition of abscissae to the 15-point gauss
c                       rule (resg).
c
c              abserr - double precison
c                       estimate of the modulus of the modulus,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void           DQK31(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief    subroutine DQK41(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 41-point
c                       gauss-kronrod rule (resk) obtained by optimal
c                       addition of abscissae to the 20-point gauss
c                       rule (resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integal of abs(f-i/(b-a))
c                       over (a,b)
*/

void          DQK41(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief      subroutine DQK51(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subroutine defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 51-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the 25-point gauss rule (resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void           DQK51(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief        subroutine DQK61(f,a,b,result,abserr,resabs,resasc)

integration rule
c        standard fortran subroutine
c        double precision version
c
c
c        parameters
c         on entry
c           f      - double precision
c                    function subprogram defining the integrand
c                    function f(x). the actual name for f needs to be
c                    declared e x t e r n a l in the calling program.
c
c           a      - double precision
c                    lower limit of integration
c
c           b      - double precision
c                    upper limit of integration
c
c         on return
c           result - double precision
c                    approximation to the integral i
c                    result is computed by applying the 61-point
c                    kronrod rule (resk) obtained by optimal addition of
c                    abscissae to the 30-point gauss rule (resg).
c
c           abserr - double precision
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed dabs(i-result)
c
c           resabs - double precision
c                    approximation to the integral j
c
c           resasc - double precision
c                    approximation to the integral of dabs(f-i/(b-a))
c
*/

void        DQK61(double(*)(double *), double *, double *, double *, double *, double *, double *);

/*
@brief      subroutine DQMOMO(alfa,beta,ri,rj,rg,rh,integr)

modified chebyshev moments
c        standard fortran subroutine
c        double precision version
c
c        parameters
c           alfa   - double precision
c                    parameter in the weight function w(x), alfa.gt.(-1)
c
c           beta   - double precision
c                    parameter in the weight function w(x), beta.gt.(-1)
c
c           ri     - double precision
c                    vector of dimension 25
c                    ri(k) is the integral over (-1,1) of
c                    (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
c
c           rj     - double precision
c                    vector of dimension 25
c                    rj(k) is the integral over (-1,1) of
c                    (1-x)**beta*t(k-1,x), k = 1, ..., 25.
c
c           rg     - double precision
c                    vector of dimension 25
c                    rg(k) is the integral over (-1,1) of
c                    (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           rh     - double precision
c                    vector of dimension 25
c                    rh(k) is the integral over (-1,1) of
c                    (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           integr - integer
c                    input parameter indicating the modified
c                    moments to be computed
c                    integr = 1 compute ri, rj
c                           = 2 compute ri, rj, rg
c                           = 3 compute ri, rj, rh
c                           = 4 compute ri, rj, rg, rh
*/

void        DQMOMO(double *, double *, double *, double *, double *, double *, int *);

/*
@brief   subroutine DQNG(f,a,b,epsabs,epsrel,result,abserr,neval,ier)

non-adaptive integration
c standard fortran subroutine
c double precision version
c
c           f      - double precision
c                    function subprogram defining the integrand function
c                    f(x). the actual name for f needs to be declared
c                    e x t e r n a l in the driver program.
c
c           a      - double precision
c                    lower limit of integration
c
c           b      - double precision
c                    upper limit of integration
c
c           epsabs - double precision
c                    absolute accuracy requested
c           epsrel - double precision
c                    relative accuracy requested
c                    if  epsabs.le.0
c                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
c                    the routine will end with ier = 6.
c
c         on return
c           result - double precision
c                    approximation to the integral i
c                    result is obtained by applying the 21-point
c                    gauss-kronrod rule (res21) obtained by optimal
c                    addition of abscissae to the 10-point gauss rule
c                    (res10), or by applying the 43-point rule (res43)
c                    obtained by optimal addition of abscissae to the
c                    21-point gauss-kronrod rule, or by applying the
c                    87-point rule (res87) obtained by optimal addition
c                    of abscissae to the 43-point rule.
c
c           abserr - double precision
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed abs(i-result)
c
c           neval  - integer
c                    number of integrand evaluations
c
c           ier    - ier = 0 normal and reliable termination of the
c                            routine. it is assumed that the requested
c                            accuracy has been achieved.
c                    ier.gt.0 abnormal termination of the routine. it is
c                            assumed that the requested accuracy has
c                            not been achieved.
c           error messages
c                    ier = 1 the maximum number of steps has been
c                            executed. the integral is probably too
c                            difficult to be calculated by dqng.
c                        = 6 the input is invalid, because
c                            epsabs.le.0 and
c                            epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
c                            result, abserr and neval are set to zero.
*/

void       DQNG(double(*)(double *), double *, double *, double *, double *, double *, double *,
	int *, int *);

/*
@brief    subroutine DQPSRT(limit,last,maxerr,ermax,elist,iord,nrmax)

ordering routine
c           standard fortran subroutine
c           double precision version
c
c           parameters (meaning at output)
c              limit  - integer
c                       maximum number of error estimates the list
c                       can contain
c
c              last   - integer
c                       number of error estimates currently in the list
c
c              maxerr - integer
c                       maxerr points to the nrmax-th largest error
c                       estimate currently in the list
c
c              ermax  - double precision
c                       nrmax-th largest error estimate
c                       ermax = elist(maxerr)
c
c              elist  - double precision
c                       vector of dimension last containing
c                       the error estimates
c
c              iord   - integer
c                       vector of dimension last, the first k elements
c                       of which contain pointers to the error
c                       estimates, such that
c                       elist(iord(1)),...,  elist(iord(k))
c                       form a decreasing sequence, with
c                       k = last if last.le.(limit/2+2), and
c                       k = limit+1-last otherwise
c
c              nrmax  - integer
c                       maxerr = iord(nrmax)
c
*/

void          DQPSRT(int *, int *, int *, double *, double *, int *, int *);

/*
@brief     subroutine QAG(f,a,b,epsabs,epsrel,key,result,abserr,neval,ier,
     *    limit,lenw,last,iwork,work)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c            f      - real
	 c                     function subprogam defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            key    - integer
	 c                     key for choice of local integration rule
	 c                     a gauss-kronrod pair is used with
	 c                       7 - 15 points if key.lt.2,
	 c                      10 - 21 points if key = 2,
	 c                      15 - 31 points if key = 3,
	 c                      20 - 41 points if key = 4,
	 c                      25 - 51 points if key = 5,
	 c                      30 - 61 points if key.gt.5.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for result and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c                      error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yield no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulaties.
	 c                             if the position of a local difficulty can
	 c                             be determined (i.e.singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.1 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set
	 c                             to zero.
	 c                             except when lenw is invalid, iwork(1),
	 c                             work(limit*2+1) and work(limit*3+1) are
	 c                             set to zero, work(1) is set to a and
	 c                             work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end with
	 c                    ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first k
	 c                    elements of which contain pointers to the error
	 c                    estimates over the subintervals, such that
	 c                    work(limit*3+iwork(1)),... , work(limit*3+iwork(k))
	 c                    form a decreasing sequence with k = last if
	 c                    last.le.(limit/2+2), and k = limit+1-last otherwise
	 c
	 c            work  - real
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left end
	 c                    points of the subintervals in the partition of
	 c                     (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain the
	 c                     right end points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last) contain
	 c                     the error estimates.
*/

void            QAG(float(*)(float *), float *, float *, float *, float *, int *, float *, float *,
	int *, int *, int *, int *, int *, int *, float *);

/*
@brief   subroutine QAGE(f,a,b,epsabs,epsrel,key,limit,result,abserr,
     *   neval,ier,alist,blist,rlist,elist,iord,last)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            key    - integer
	 c                     key for choice of local integration rule
	 c                     a gauss-kronrod pair is used with
	 c                          7 - 15 points if key.lt.2,
	 c                         10 - 21 points if key = 2,
	 c                         15 - 31 points if key = 3,
	 c                         20 - 41 points if key = 4,
	 c                         25 - 51 points if key = 5,
	 c                         30 - 61 points if key.gt.5.
	 c
	 c            limit  - integer
	 c                     gives an upperbound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for result and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value
	 c                             of limit.
	 c                             however, if this yields no improvement it
	 c                             is rather advised to analyze the integrand
	 c                             in order to determine the integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined(e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling the integrator on the
	 c                             subranges. if possible, an appropriate
	 c                             special-purpose integrator should be used
	 c                             which is designed for handling the type of
	 c                             difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             result, abserr, neval, last, rlist(1) ,
	 c                             elist(1) and iord(1) are set to zero.
	 c                             alist(1) and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the left
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            blist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the right
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            rlist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the
	 c                      integral approximations on the subintervals
	 c
	 c            elist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the moduli of the
	 c                      absolute error estimates on the subintervals
	 c
	 c            iord    - integer
	 c                      vector of dimension at least limit, the first k
	 c                      elements of which are pointers to the
	 c                      error estimates over the subintervals,
	 c                      such that elist(iord(1)), ...,
	 c                      elist(iord(k)) form a decreasing sequence,
	 c                      with k = last if last.le.(limit/2+2), and
	 c                      k = limit+1-last otherwise
	 c
	 c            last    - integer
	 c                      number of subintervals actually produced in the
	 c                      subdivision process
*/

void             QAGE(float(*)(float *), float *, float *, float *, float *, int *, int *, float *, float *,
	int *, int *, float *, float *, float *, float *, int *, int *);

/*
@brief  subroutine QAGI(f,bound,inf,epsabs,epsrel,result,abserr,neval,
     *   ier,limit,lenw,last,iwork,work)

	 integration over infinite intervals
	 c        standard fortran subroutine
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            bound  - real
	 c                     finite bound of integration range
	 c                     (has no meaning if interval is doubly-infinite)
	 c
	 c            inf    - integer
	 c                     indicating the kind of integration range involved
	 c                     inf = 1 corresponds to  (bound,+infinity),
	 c                     inf = -1            to  (-infinity,bound),
	 c                     inf = 2             to (-infinity,+infinity).
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for result and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is assumed that the requested tolerance
	 c                             cannot be achieved, and that the returned
	 c                             result is the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                              or limit.lt.1 or leniw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when limit or leniw is
	 c                             invalid, iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first
	 c                    k elements of which contain pointers
	 c                    to the error estimates over the subintervals,
	 c                    such that work(limit*3+iwork(1)),... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2), and
	 c                    k = limit+1-last otherwise
	 c
	 c            work  - real
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ...,work(limit*2+last) contain the
	 c                     integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3)
	 c                     contain the error estimates.
*/

void            QAGI(float(*)(float *), float *, int *, float *, float *, float *, float *, int *, int *,
	int *, int *, int *, int *, float *);

/*
@brief  subroutine QAGIE(f,bound,inf,epsabs,epsrel,limit,result,abserr,
     *   neval,ier,alist,blist,rlist,elist,iord,last)

	 integration over infinite intervals
	 c standard fortran subroutine
	 c
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            bound  - real
	 c                     finite bound of integration range
	 c                     (has no meaning if interval is doubly-infinite)
	 c
	 c            inf    - real
	 c                     indicating the kind of integration range involved
	 c                     inf = 1 corresponds to  (bound,+infinity),
	 c                     inf = -1            to  (-infinity,bound),
	 c                     inf = 2             to (-infinity,+infinity).
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for result and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however,if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty can
	 c                             be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is assumed that the requested tolerance
	 c                             cannot be achieved, and that the returned
	 c                             result is the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             result, abserr, neval, last, rlist(1),
	 c                             elist(1) and iord(1) are set to zero.
	 c                             alist(1) and blist(1) are set to 0
	 c                             and 1 respectively.
	 c
	 c            alist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the transformed integration range (0,1).
	 c
	 c            blist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the transformed integration range (0,1).
	 c
	 c            rlist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - real
	 c                     vector of dimension at least limit,  the first
	 c                     last elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced
	 c                     in the subdivision process
*/

void             QAGIE(float(*)(float *), float *, float *, float *, float *, int *, float *, float *,
	int *, int *, float *, float *, float *, float *, int *, int *);

/*
@brief       subroutine QAGP(f,a,b,npts2,points,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,lenw,last,iwork,work)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            npts2  - integer
	 c                     number equal to two more than the number of
	 c                     user-supplied break points within the integration
	 c                     range, npts.ge.2.
	 c                     if npts2.lt.2, the routine will end with ier = 6.
	 c
	 c            points - real
	 c                     vector of dimension npts2, the first (npts2-2)
	 c                     elements of which are the user provided break
	 c                     points. if these points do not constitute an
	 c                     ascending sequence there will be an automatic
	 c                     sorting.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (i.e. singularity,
	 c                             discontinuity within the interval), it
	 c                             should be supplied to the routine as an
	 c                             element of the vector points. if necessary
	 c                             an appropriate special-purpose integrator
	 c                             must be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved, and that
	 c                             the returned result is the best which
	 c                             can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid because
	 c                             npts2.lt.2 or
	 c                             break points are specified outside
	 c                             the integration range or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when leniw or lenw or npts2 is
	 c                             invalid, iwork(1), iwork(limit+1),
	 c                             work(limit*2+1) and work(limit*3+1)
	 c                             are set to zero.
	 c                             work(1) is set to a and work(limit+1)
	 c                             to b (where limit = (leniw-npts2)/2).
	 c
	 c         dimensioning parameters
	 c            leniw - integer
	 c                    dimensioning parameter for iwork
	 c                    leniw determines limit = (leniw-npts2)/2,
	 c                    which is the maximum number of subintervals in the
	 c                    partition of the given integration interval (a,b),
	 c                    leniw.ge.(3*npts2-2).
	 c                    if leniw.lt.(3*npts2-2), the routine will end with
	 c                    ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least leniw*2-npts2.
	 c                    if lenw.lt.leniw*2-npts2, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least leniw. on return,
	 c                    the first k elements of which contain
	 c                    pointers to the error estimates over the
	 c                    subintervals, such that work(limit*3+iwork(1)),...,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2), and
	 c                    k = limit+1-last otherwise
	 c                    iwork(limit+1), ...,iwork(limit+last) contain the
	 c                     subdivision levels of the subintervals, i.e.
	 c                     if (aa,bb) is a subinterval of (p1,p2)
	 c                     where p1 as well as p2 is a user-provided
	 c                     break point or integration limit, then (aa,bb) has
	 c                     level l if abs(bb-aa) = abs(p2-p1)*2**(-l),
	 c                    iwork(limit*2+1), ..., iwork(limit*2+npts2) have
	 c                     no significance for the user,
	 c                    note that limit = (leniw-npts2)/2.
	 c
	 c            work  - real
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the corresponding error estimates,
	 c                    work(limit*4+1), ..., work(limit*4+npts2)
	 c                     contain the integration limits and the
	 c                     break points sorted in an ascending sequence.
	 c                    note that limit = (leniw-npts2)/2.
*/

void                 QAGP(float(*)(float *), float *, float *, int *, float *, float *, float *, float *,
	float *, int *, int *, int *, int *, int *, int *, float *);

/*
@brief   subroutine QAGPE(f,a,b,npts2,points,epsabs,epsrel,limit,result,
     *   abserr,neval,ier,alist,blist,rlist,elist,pts,iord,level,ndin,
     *   last)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            npts2  - integer
	 c                     number equal to two more than the number of
	 c                     user-supplied break points within the integration
	 c                     range, npts2.ge.2.
	 c                     if npts2.lt.2, the routine will end with ier = 6.
	 c
	 c            points - real
	 c                     vector of dimension npts2, the first (npts2-2)
	 c                     elements of which are the user provided break
	 c                     points. if these points do not constitute an
	 c                     ascending sequence there will be an automatic
	 c                     sorting.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.npts2
	 c                     if limit.lt.npts2, the routine will end with
	 c                     ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (i.e. singularity,
	 c                             discontinuity within the interval), it
	 c                             should be supplied to the routine as an
	 c                             element of the vector points. if necessary
	 c                             an appropriate special-purpose integrator
	 c                             must be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be
	 c                             achieved, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid because
	 c                             npts2.lt.2 or
	 c                             break points are specified outside
	 c                             the integration range or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.npts2.
	 c                             result, abserr, neval, last, rlist(1),
	 c                             and elist(1) are set to zero. alist(1) and
	 c                             blist(1) are set to a and b respectively.
	 c
	 c            alist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            blist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            rlist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            pts    - real
	 c                     vector of dimension at least npts2, containing the
	 c                     integration limits and the break points of the
	 c                     interval in ascending sequence.
	 c
	 c            level  - integer
	 c                     vector of dimension at least limit, containing the
	 c                     subdivision levels of the subinterval, i.e. if
	 c                     (aa,bb) is a subinterval of (p1,p2) where p1 as
	 c                     well as p2 is a user-provided break point or
	 c                     integration limit, then (aa,bb) has level l if
	 c                     abs(bb-aa) = abs(p2-p1)*2**(-l).
	 c
	 c            ndin   - integer
	 c                     vector of dimension at least npts2, after first
	 c                     integration over the intervals (pts(i)),pts(i+1),
	 c                     i = 0,1, ..., npts2-2, the error estimates over
	 c                     some of the intervals may have been increased
	 c                     artificially, in order to put their subdivision
	 c                     forward. if this happens for the subinterval
	 c                     numbered k, ndin(k) is put to 1, otherwise
	 c                     ndin(k) = 0.
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in the
	 c                     subdivisions process
*/

void             QAGPE(float(*)(float *), float *, float *, int *, float *, float *, float *, int *, float *,
	float *, int *, int *, float *, float *, float *, float *, float *, int *,
	int *, int *, int *);

/*
@brief    subroutine QAGS(f,a,b,epsabs,epsrel,result,abserr,neval,ier,
     *   limit,lenw,last,iwork,work)

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account. however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be
	 c                             achieved, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28)
	 c                             or limit.lt.1 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero.except when limit or lenw is invalid,
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw  - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end
	 c                    with ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, detemines the
	 c                    number of significant elements actually in the work
	 c                    arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first k
	 c                    elements of which contain pointers
	 c                    to the error estimates over the subintervals
	 c                    such that work(limit*3+iwork(1)),... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2),
	 c                    and k = limit+1-last otherwise
	 c
	 c            work  - real
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end-points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end-points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the error estimates.
	 c
*/

void              QAGS(float(*)(float *), float *, float *, float *, float *, float *, float *,
	int *, int *, int *, int *, int *, int *, float *);

/*
@brief      subroutine QAGSE(f,a,b,epsabs,epsrel,limit,result,abserr,neval,
     *   ier,alist,blist,rlist,elist,iord,last)   

	 computation of a definite integral
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upperbound on the number of subintervals
	 c                     in the partition of (a,b)
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                         = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties. if
	 c                             the position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used, which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved, and that the
	 c                             returned result is the best which can be
	 c                             obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             epsabs.le.0 and
	 c                             epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
	 c                             result, abserr, neval, last, rlist(1),
	 c                             iord(1) and elist(1) are set to zero.
	 c                             alist(1) and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left end points
	 c                     of the subintervals in the partition of the
	 c                     given integration range (a,b)
	 c
	 c            blist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right end points
	 c                     of the subintervals in the partition of the given
	 c                     integration range (a,b)
	 c
	 c            rlist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the
	 c                     error estimates over the subintervals,
	 c                     such that elist(iord(1)), ..., elist(iord(k))
	 c                     form a decreasing sequence, with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in the
	 c                     subdivision process
*/

void              QAGSE(float(*)(float *), float *, float *, float *, float *, int *, float *, float *,
	int *, int *, float *, float *, float *, float *, int *, int *);

/*
@brief    subroutine QAWC(f,a,b,c,epsabs,epsrel,result,abserr,neval,ier,
     *   limit,lenw,last,iwork,work)

	 computation of a cauchy principal value
	 c        standard fortran subroutine
	 c        real version
	 c
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     under limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            c      - parameter in the weight function, c.ne.a, c.ne.b.
	 c                     if c = a or c = b, the routine will end with
	 c                     ier = 6 .
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate or the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of limit
	 c                             (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty
	 c                             can be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling
	 c                             appropriate integrators on the subranges.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             c = a or c = b or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.1 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. exept when lenw or limit is invalid,
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit - integer
	 c                    dimensioning parameter for iwork
	 c                    limit determines the maximum number of subintervals
	 c                    in the partition of the given integration interval
	 c                    (a,b), limit.ge.1.
	 c                    if limit.lt.1, the routine will end with ier = 6.
	 c
	 c           lenw   - integer
	 c                    dimensioning parameter for work
	 c                    lenw must be at least limit*4.
	 c                    if lenw.lt.limit*4, the routine will end with
	 c                    ier = 6.
	 c
	 c            last  - integer
	 c                    on return, last equals the number of subintervals
	 c                    produced in the subdivision process, which
	 c                    determines the number of significant elements
	 c                    actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork - integer
	 c                    vector of dimension at least limit, the first k
	 c                    elements of which contain pointers
	 c                    to the error estimates over the subintervals,
	 c                    such that work(limit*3+iwork(1)), ... ,
	 c                    work(limit*3+iwork(k)) form a decreasing
	 c                    sequence, with k = last if last.le.(limit/2+2),
	 c                    and k = limit+1-last otherwise
	 c
	 c            work  - real
	 c                    vector of dimension at least lenw
	 c                    on return
	 c                    work(1), ..., work(last) contain the left
	 c                     end points of the subintervals in the
	 c                     partition of (a,b),
	 c                    work(limit+1), ..., work(limit+last) contain
	 c                     the right end points,
	 c                    work(limit*2+1), ..., work(limit*2+last) contain
	 c                     the integral approximations over the subintervals,
	 c                    work(limit*3+1), ..., work(limit*3+last)
	 c                     contain the error estimates.
*/

void             QAWC(float(*)(float *), float *, float *, float *, float *, float *, float *,
	float *, int *, int *, int *, int *, int *, int *, float *);

/*
@brief    subroutine QAWCE(f,a,b,c,epsabs,epsrel,limit,result,abserr,neval,
     *   ier,alist,blist,rlist,elist,iord,last)

	 computation of a cauchy principal value
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            c      - real
	 c                     parameter in the weight function, c.ne.a, c.ne.b
	 c                     if c = a or c = b, the routine will end with
	 c                     ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.1
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more sub-
	 c                             divisions by increasing the value of
	 c                             limit. however, if this yields no
	 c                             improvement it is advised to analyze the
	 c                             the integrand, in order to determine the
	 c                             the integration difficulties. if the
	 c                             position of a local difficulty can be
	 c                             determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling
	 c                             appropriate integrators on the subranges.
	 c                         = 2 the occurrence of roundoff error is detec-
	 c                             ted, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour
	 c                             occurs at some interior points of
	 c                             the integration interval.
	 c                         = 6 the input is invalid, because
	 c                             c = a or c = b or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.1.
	 c                             result, abserr, neval, rlist(1), elist(1),
	 c                             iord(1) and last are set to zero. alist(1)
	 c                             and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the left
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            blist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the right
	 c                      end points of the subintervals in the partition
	 c                      of the given integration range (a,b)
	 c
	 c            rlist   - real
	 c                      vector of dimension at least limit, the first
	 c                       last  elements of which are the integral
	 c                      approximations on the subintervals
	 c
	 c            elist   - real
	 c                      vector of dimension limit, the first  last
	 c                      elements of which are the moduli of the absolute
	 c                      error estimates on the subintervals
	 c
	 c            iord    - integer
	 c                      vector of dimension at least limit, the first k
	 c                      elements of which are pointers to the error
	 c                      estimates over the subintervals, so that
	 c                      elist(iord(1)), ..., elist(iord(k)) with k = last
	 c                      if last.le.(limit/2+2), and k = limit+1-last
	 c                      otherwise, form a decreasing sequence
	 c
	 c            last    - integer
	 c                      number of subintervals actually produced in
	 c                      the subdivision process
*/

void            QAWCE(float(*)(float *), float *, float *, float *, float *, float *, int *,
	float *, float *, int *, int *, float *, float *, float *, float *, int *, int *);

/*
@brief        subroutine QAWFE(f,a,omega,integr,epsabs,limlst,limit,maxp1,
     *   result,abserr,neval,ier,rslst,erlst,ierlst,lst,alist,blist,
     *   rlist,elist,iord,nnlog,chebmo)

	 computation of fourier integrals
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to
	 c                     be declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            omega  - real
	 c                     parameter in the weight function
	 c
	 c            integr - integer
	 c                     indicates which weight function is used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1.and.integr.ne.2, the routine will
	 c                     end with ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested, epsabs.gt.0
	 c                     if epsabs.le.0, the routine will end with ier = 6.
	 c
	 c            limlst - integer
	 c                     limlst gives an upper bound on the number of
	 c                     cycles, limlst.ge.1.
	 c                     if limlst.lt.3, the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     allowed in the partition of each cycle, limit.ge.1
	 c                     each cycle, limit.ge.1.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of
	 c                     chebyshev moments which can be stored, i.e.
	 c                     for the intervals of lengths abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral x
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - ier = 0 normal and reliable termination of
	 c                             the routine. it is assumed that the
	 c                             requested accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine. the
	 c                             estimates for integral and error are less
	 c                             reliable. it is assumed that the requested
	 c                             accuracy has not been achieved.
	 c            error messages
	 c                    if omega.ne.0
	 c                     ier = 1 maximum number of  cycles  allowed
	 c                             has been achieved., i.e. of subintervals
	 c                             (a+(k-1)c,a+kc) where
	 c                             c = (2*int(abs(omega))+1)*pi/abs(omega),
	 c                             for k = 1, 2, ..., lst.
	 c                             one can allow more cycles by increasing
	 c                             the value of limlst (and taking the
	 c                             according dimension adjustments into
	 c                             account).
	 c                             examine the array iwork which contains
	 c                             the error flags on the cycles, in order to
	 c                             look for eventual local integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined (e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling appropriate integrators on
	 c                             the subranges.
	 c                         = 4 the extrapolation table constructed for
	 c                             convergence acceleration of the series
	 c                             formed by the integral contributions over
	 c                             the cycles, does not converge to within
	 c                             the requested accuracy. as in the case of
	 c                             ier = 1, it is advised to examine the
	 c                             array iwork which contains the error
	 c                             flags on the cycles.
	 c                         = 6 the input is invalid because
	 c                             (integr.ne.1 and integr.ne.2) or
	 c                              epsabs.le.0 or limlst.lt.3.
	 c                              result, abserr, neval, lst are set
	 c                              to zero.
	 c                         = 7 bad integrand behaviour occurs within one
	 c                             or more of the cycles. location and type
	 c                             of the difficulty involved can be
	 c                             determined from the vector ierlst. here
	 c                             lst is the number of cycles actually
	 c                             needed (see below).
	 c                             ierlst(k) = 1 the maximum number of
	 c                                           subdivisions (= limit) has
	 c                                           been achieved on the k th
	 c                                           cycle.
	 c                                       = 2 occurrence of roundoff error
	 c                                           is detected and prevents the
	 c                                           tolerance imposed on the
	 c                                           k th cycle, from being
	 c                                           achieved.
	 c                                       = 3 extremely bad integrand
	 c                                           behaviour occurs at some
	 c                                           points of the k th cycle.
	 c                                       = 4 the integration procedure
	 c                                           over the k th cycle does
	 c                                           not converge (to within the
	 c                                           required accuracy) due to
	 c                                           roundoff in the
	 c                                           extrapolation procedure
	 c                                           invoked on this cycle. it
	 c                                           is assumed that the result
	 c                                           on this interval is the
	 c                                           best which can be obtained.
	 c                                       = 5 the integral over the k th
	 c                                           cycle is probably divergent
	 c                                           or slowly convergent. it
	 c                                           must be noted that
	 c                                           divergence can occur with
	 c                                           any other value of
	 c                                           ierlst(k).
	 c                    if omega = 0 and integr = 1,
	 c                    the integral is calculated by means of dqagie
	 c                    and ier = ierlst(1) (with meaning as described
	 c                    for ierlst(k), k = 1).
	 c
	 c            rslst  - real
	 c                     vector of dimension at least limlst
	 c                     rslst(k) contains the integral contribution
	 c                     over the interval (a+(k-1)c,a+kc) where
	 c                     c = (2*int(abs(omega))+1)*pi/abs(omega),
	 c                     k = 1, 2, ..., lst.
	 c                     note that, if omega = 0, rslst(1) contains
	 c                     the value of the integral over (a,infinity).
	 c
	 c            erlst  - real
	 c                     vector of dimension at least limlst
	 c                     erlst(k) contains the error estimate corresponding
	 c                     with rslst(k).
	 c
	 c            ierlst - integer
	 c                     vector of dimension at least limlst
	 c                     ierlst(k) contains the error flag corresponding
	 c                     with rslst(k). for the meaning of the local error
	 c                     flags see description of output parameter ier.
	 c
	 c            lst    - integer
	 c                     number of subintervals needed for the integration
	 c                     if omega = 0 then lst is set to 1.
	 c
	 c            alist, blist, rlist, elist - real
	 c                     vector of dimension at least limit,
	 c
	 c            iord, nnlog - integer
	 c                     vector of dimension at least limit, providing
	 c                     space for the quantities needed in the subdivision
	 c                     process of each cycle
	 c
	 c            chebmo - real
	 c                     array of dimension at least (maxp1,25), providing
	 c                     space for the chebyshev moments needed within the
	 c                     cycles
*/

void             QAWFE(float(*)(float *), float *, float *, int *, float *, int *, int *, int *,
	float *, float *, int *, int *, float *, float *, int *, int *, float *, float *, float *,
	float *, int *, int *, float *);

/*
@brief  subroutine QAWO(f,a,b,omega,integr,epsabs,epsrel,result,abserr,
     *   neval,ier,leniw,maxp1,lenw,last,iwork,work)

	 computation of oscillatory integrals
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the function
	 c                     f(x).  the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            omega  - real
	 c                     parameter in the integrand weight function
	 c
	 c            integr - integer
	 c                     indicates which of the weight functions is used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1.and.integr.ne.2, the routine will
	 c                     end with ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if epsabs.le.0 and
	 c                     epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of  integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             (= leniw/2) has been achieved. one can
	 c                             allow more subdivisions by increasing the
	 c                             value of leniw (and taking the according
	 c                             dimension adjustments into account).
	 c                             however, if this yields no improvement it
	 c                             is advised to analyze the integrand in
	 c                             order to determine the integration
	 c                             difficulties. if the position of a local
	 c                             difficulty can be determined (e.g.
	 c                             singularity, discontinuity within the
	 c                             interval) one will probably gain from
	 c                             splitting up the interval at this point
	 c                             and calling the integrator on the
	 c                             subranges. if possible, an appropriate
	 c                             special-purpose integrator should be used
	 c                             which is designed for handling the type of
	 c                             difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some interior points of the
	 c                             integration interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table. it is presumed that
	 c                             the requested tolerance cannot be achieved
	 c                             due to roundoff in the extrapolation
	 c                             table, and that the returned result is
	 c                             the best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or (integr.ne.1 and integr.ne.2),
	 c                             or leniw.lt.2 or maxp1.lt.1 or
	 c                             lenw.lt.leniw*2+maxp1*25.
	 c                             result, abserr, neval, last are set to
	 c                             zero. except when leniw, maxp1 or lenw are
	 c                             invalid, work(limit*2+1), work(limit*3+1),
	 c                             iwork(1), iwork(limit+1) are set to zero,
	 c                             work(1) is set to a and work(limit+1) to
	 c                             b.
	 c
	 c         dimensioning parameters
	 c            leniw  - integer
	 c                     dimensioning parameter for iwork.
	 c                     leniw/2 equals the maximum number of subintervals
	 c                     allowed in the partition of the given integration
	 c                     interval (a,b), leniw.ge.2.
	 c                     if leniw.lt.2, the routine will end with ier = 6.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of chebyshev
	 c                     moments which can be stored, i.e. for the
	 c                     intervals of lengths abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1
	 c                     if maxp1.lt.1, the routine will end with ier = 6.
	 c
	 c            lenw   - integer
	 c                     dimensioning parameter for work
	 c                     lenw must be at least leniw*2+maxp1*25.
	 c                     if lenw.lt.(leniw*2+maxp1*25), the routine will
	 c                     end with ier = 6.
	 c
	 c            last   - integer
	 c                     on return, last equals the number of subintervals
	 c                     produced in the subdivision process, which
	 c                     determines the number of significant elements
	 c                     actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork  - integer
	 c                     vector of dimension at least leniw
	 c                     on return, the first k elements of which contain
	 c                     pointers to the error estimates over the
	 c                     subintervals, such that work(limit*3+iwork(1)), ..
	 c                     work(limit*3+iwork(k)) form a decreasing
	 c                     sequence, with limit = lenw/2 , and k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise.
	 c                     furthermore, iwork(limit+1), ..., iwork(limit+
	 c                     last) indicate the subdivision levels of the
	 c                     subintervals, such that iwork(limit+i) = l means
	 c                     that the subinterval numbered i is of length
	 c                     abs(b-a)*2**(1-l).
	 c
	 c            work   - real
	 c                     vector of dimension at least lenw
	 c                     on return
	 c                     work(1), ..., work(last) contain the left
	 c                      end points of the subintervals in the
	 c                      partition of (a,b),
	 c                     work(limit+1), ..., work(limit+last) contain
	 c                      the right end points,
	 c                     work(limit*2+1), ..., work(limit*2+last) contain
	 c                      the integral approximations over the
	 c                      subintervals,
	 c                     work(limit*3+1), ..., work(limit*3+last)
	 c                      contain the error estimates.
	 c                     work(limit*4+1), ..., work(limit*4+maxp1*25)
	 c                      provide space for storing the chebyshev moments.
	 c                     note that limit = lenw/2.
*/


void             QAWO(float(*)(float *), float *, float *, float *, int *, float *, float *,
	float *, float *, int *, int *, int *, int *, int *, int *, int *, float *);

/*
@brief   subroutine  QAWOE(f,a,b,omega,integr,epsabs,epsrel,limit,icall,
     *  maxp1,result,abserr,neval,ier,last,alist,blist,rlist,elist,iord,
     *   nnlog,momcom,chebmo)

	 computation of oscillatory integrals
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration
	 c
	 c            omega  - real
	 c                     parameter in the integrand weight function
	 c
	 c            integr - integer
	 c                     indicates which of the weight functions is to be
	 c                     used
	 c                     integr = 1      w(x) = cos(omega*x)
	 c                     integr = 2      w(x) = sin(omega*x)
	 c                     if integr.ne.1 and integr.ne.2, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subdivisions
	 c                     in the partition of (a,b), limit.ge.1.
	 c
	 c            icall  - integer
	 c                     if dqawoe is to be used only once, icall must
	 c                     be set to 1.  assume that during this call, the
	 c                     chebyshev moments (for clenshaw-curtis integration
	 c                     of degree 24) have been computed for intervals of
	 c                     lenghts (abs(b-a))*2**(-l), l=0,1,2,...momcom-1.
	 c                     if icall.gt.1 this means that dqawoe has been
	 c                     called twice or more on intervals of the same
	 c                     length abs(b-a). the chebyshev moments already
	 c                     computed are then re-used in subsequent calls.
	 c                     if icall.lt.1, the routine will end with ier = 6.
	 c
	 c            maxp1  - integer
	 c                     gives an upper bound on the number of chebyshev
	 c                     moments which can be stored, i.e. for the
	 c                     intervals of lenghts abs(b-a)*2**(-l),
	 c                     l=0,1, ..., maxp1-2, maxp1.ge.1.
	 c                     if maxp1.lt.1, the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the
	 c                             requested accuracy has been achieved.
	 c                   - ier.gt.0 abnormal termination of the routine.
	 c                             the estimates for integral and error are
	 c                             less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand, in order to
	 c                             determine the integration difficulties.
	 c                             if the position of a local difficulty can
	 c                             be determined (e.g. singularity,
	 c                             discontinuity within the interval) one
	 c                             will probably gain from splitting up the
	 c                             interval at this point and calling the
	 c                             integrator on the subranges. if possible,
	 c                             an appropriate special-purpose integrator
	 c                             should be used which is designed for
	 c                             handling the type of difficulty involved.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                             the error may be under-estimated.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 4 the algorithm does not converge.
	 c                             roundoff error is detected in the
	 c                             extrapolation table.
	 c                             it is presumed that the requested
	 c                             tolerance cannot be achieved due to
	 c                             roundoff in the extrapolation table,
	 c                             and that the returned result is the
	 c                             best which can be obtained.
	 c                         = 5 the integral is probably divergent, or
	 c                             slowly convergent. it must be noted that
	 c                             divergence can occur with any other value
	 c                             of ier.gt.0.
	 c                         = 6 the input is invalid, because
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or (integr.ne.1 and integr.ne.2) or
	 c                             icall.lt.1 or maxp1.lt.1.
	 c                             result, abserr, neval, last, rlist(1),
	 c                             elist(1), iord(1) and nnlog(1) are set
	 c                             to zero. alist(1) and blist(1) are set
	 c                             to a and b respectively.
	 c
	 c            last  -  integer
	 c                     on return, last equals the number of
	 c                     subintervals produces in the subdivision
	 c                     process, which determines the number of
	 c                     significant elements actually in the
	 c                     work arrays.
	 c            alist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            blist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            rlist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     elements of which are pointers to the error
	 c                     estimates over the subintervals,
	 c                     such that elist(iord(1)), ...,
	 c                     elist(iord(k)) form a decreasing sequence, with
	 c                     k = last if last.le.(limit/2+2), and
	 c                     k = limit+1-last otherwise.
	 c
	 c            nnlog  - integer
	 c                     vector of dimension at least limit, containing the
	 c                     subdivision levels of the subintervals, i.e.
	 c                     iwork(i) = l means that the subinterval
	 c                     numbered i is of length abs(b-a)*2**(1-l)
	 c
	 c         on entry and return
	 c            momcom - integer
	 c                     indicating that the chebyshev moments
	 c                     have been computed for intervals of lengths
	 c                     (abs(b-a))*2**(-l), l=0,1,2, ..., momcom-1,
	 c                     momcom.lt.maxp1
	 c
	 c            chebmo - real
	 c                     array of dimension (maxp1,25) containing the
	 c                     chebyshev moments
*/

void            QAWOE(float(*)(float *), float *, float *, float *, int *, float *, float *, int *, int *,
	int *, float *, float *, int *, int *, int *, float *, float *, float *, float *,
	int *, int *, int *, float *);

/*
@brief  subroutine QAWS(f,a,b,alfa,beta,integr,epsabs,epsrel,result,
     *   abserr,neval,ier,limit,lenw,last,iwork,work)

	 integration of functions having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration, b.gt.a
	 c                     if b.le.a, the routine will end with ier = 6.
	 c
	 c            alfa   - real
	 c                     parameter in the integrand function, alfa.gt.(-1)
	 c                     if alfa.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            beta   - real
	 c                     parameter in the integrand function, beta.gt.(-1)
	 c                     if beta.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            integr - integer
	 c                     indicates which weight function is to be used
	 c                     = 1  (x-a)**alfa*(b-x)**beta
	 c                     = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                     = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                     = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
	 c                     if integr.lt.1 or integr.gt.4, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for the integral and error
	 c                             are less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                     ier = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit (and taking the according dimension
	 c                             adjustments into account). however, if
	 c                             this yields no improvement it is advised
	 c                             to analyze the integrand, in order to
	 c                             determine the integration difficulties
	 c                             which prevent the requested tolerance from
	 c                             being achieved. in case of a jump
	 c                             discontinuity or a local singularity
	 c                             of algebraico-logarithmic type at one or
	 c                             more interior points of the integration
	 c                             range, one should proceed by splitting up
	 c                             the interval at these points and calling
	 c                             the integrator on the subranges.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             b.le.a or alfa.le.(-1) or beta.le.(-1) or
	 c                             or integr.lt.1 or integr.gt.4 or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
	 c                             or limit.lt.2 or lenw.lt.limit*4.
	 c                             result, abserr, neval, last are set to
	 c                             zero. except when lenw or limit is invalid
	 c                             iwork(1), work(limit*2+1) and
	 c                             work(limit*3+1) are set to zero, work(1)
	 c                             is set to a and work(limit+1) to b.
	 c
	 c         dimensioning parameters
	 c            limit  - integer
	 c                     dimensioning parameter for iwork
	 c                     limit determines the maximum number of
	 c                     subintervals in the partition of the given
	 c                     integration interval (a,b), limit.ge.2.
	 c                     if limit.lt.2, the routine will end with ier = 6.
	 c
	 c            lenw   - integer
	 c                     dimensioning parameter for work
	 c                     lenw must be at least limit*4.
	 c                     if lenw.lt.limit*4, the routine will end
	 c                     with ier = 6.
	 c
	 c            last   - integer
	 c                     on return, last equals the number of
	 c                     subintervals produced in the subdivision process,
	 c                     which determines the significant number of
	 c                     elements actually in the work arrays.
	 c
	 c         work arrays
	 c            iwork  - integer
	 c                     vector of dimension limit, the first k
	 c                     elements of which contain pointers
	 c                     to the error estimates over the subintervals,
	 c                     such that work(limit*3+iwork(1)), ...,
	 c                     work(limit*3+iwork(k)) form a decreasing
	 c                     sequence with k = last if last.le.(limit/2+2),
	 c                     and k = limit+1-last otherwise
	 c
	 c            work   - real
	 c                     vector of dimension lenw
	 c                     on return
	 c                     work(1), ..., work(last) contain the left
	 c                      end points of the subintervals in the
	 c                      partition of (a,b),
	 c                     work(limit+1), ..., work(limit+last) contain
	 c                      the right end points,
	 c                     work(limit*2+1), ..., work(limit*2+last)
	 c                      contain the integral approximations over
	 c                      the subintervals,
	 c                     work(limit*3+1), ..., work(limit*3+last)
	 c                      contain the error estimates.
*/

void            QAWS(float(*)(float *), float *, float *, float *, float *, int *, float *, float *,
	float *, float *, int *, int *, int *, int *, int *, int *, float *);

/*
@brief    subroutine QAWSE(f,a,b,alfa,beta,integr,epsabs,epsrel,limit,
     *   result,abserr,neval,ier,alist,blist,rlist,elist,iord,last)

	 integration of functions having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c         on entry
	 c            f      - real
	 c                     function subprogram defining the integrand
	 c                     function f(x). the actual name for f needs to be
	 c                     declared e x t e r n a l in the driver program.
	 c
	 c            a      - real
	 c                     lower limit of integration
	 c
	 c            b      - real
	 c                     upper limit of integration, b.gt.a
	 c                     if b.le.a, the routine will end with ier = 6.
	 c
	 c            alfa   - real
	 c                     parameter in the weight function, alfa.gt.(-1)
	 c                     if alfa.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            beta   - real
	 c                     parameter in the weight function, beta.gt.(-1)
	 c                     if beta.le.(-1), the routine will end with
	 c                     ier = 6.
	 c
	 c            integr - integer
	 c                     indicates which weight function is to be used
	 c                     = 1  (x-a)**alfa*(b-x)**beta
	 c                     = 2  (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                     = 3  (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                     = 4  (x-a)**alfa*(b-x)**beta*log(x-a)*log(b-x)
	 c                     if integr.lt.1 or integr.gt.4, the routine
	 c                     will end with ier = 6.
	 c
	 c            epsabs - real
	 c                     absolute accuracy requested
	 c            epsrel - real
	 c                     relative accuracy requested
	 c                     if  epsabs.le.0
	 c                     and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                     the routine will end with ier = 6.
	 c
	 c            limit  - integer
	 c                     gives an upper bound on the number of subintervals
	 c                     in the partition of (a,b), limit.ge.2
	 c                     if limit.lt.2, the routine will end with ier = 6.
	 c
	 c         on return
	 c            result - real
	 c                     approximation to the integral
	 c
	 c            abserr - real
	 c                     estimate of the modulus of the absolute error,
	 c                     which should equal or exceed abs(i-result)
	 c
	 c            neval  - integer
	 c                     number of integrand evaluations
	 c
	 c            ier    - integer
	 c                     ier = 0 normal and reliable termination of the
	 c                             routine. it is assumed that the requested
	 c                             accuracy has been achieved.
	 c                     ier.gt.0 abnormal termination of the routine
	 c                             the estimates for the integral and error
	 c                             are less reliable. it is assumed that the
	 c                             requested accuracy has not been achieved.
	 c            error messages
	 c                         = 1 maximum number of subdivisions allowed
	 c                             has been achieved. one can allow more
	 c                             subdivisions by increasing the value of
	 c                             limit. however, if this yields no
	 c                             improvement, it is advised to analyze the
	 c                             integrand in order to determine the
	 c                             integration difficulties which prevent the
	 c                             requested tolerance from being achieved.
	 c                             in case of a jump discontinuity or a local
	 c                             singularity of algebraico-logarithmic type
	 c                             at one or more interior points of the
	 c                             integration range, one should proceed by
	 c                             splitting up the interval at these
	 c                             points and calling the integrator on the
	 c                             subranges.
	 c                         = 2 the occurrence of roundoff error is
	 c                             detected, which prevents the requested
	 c                             tolerance from being achieved.
	 c                         = 3 extremely bad integrand behaviour occurs
	 c                             at some points of the integration
	 c                             interval.
	 c                         = 6 the input is invalid, because
	 c                             b.le.a or alfa.le.(-1) or beta.le.(-1), or
	 c                             integr.lt.1 or integr.gt.4, or
	 c                             (epsabs.le.0 and
	 c                              epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
	 c                             or limit.lt.2.
	 c                             result, abserr, neval, rlist(1), elist(1),
	 c                             iord(1) and last are set to zero. alist(1)
	 c                             and blist(1) are set to a and b
	 c                             respectively.
	 c
	 c            alist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the left
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            blist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the right
	 c                     end points of the subintervals in the partition
	 c                     of the given integration range (a,b)
	 c
	 c            rlist  - real
	 c                     vector of dimension at least limit,the first
	 c                      last  elements of which are the integral
	 c                     approximations on the subintervals
	 c
	 c            elist  - real
	 c                     vector of dimension at least limit, the first
	 c                      last  elements of which are the moduli of the
	 c                     absolute error estimates on the subintervals
	 c
	 c            iord   - integer
	 c                     vector of dimension at least limit, the first k
	 c                     of which are pointers to the error
	 c                     estimates over the subintervals, so that
	 c                     elist(iord(1)), ..., elist(iord(k)) with k = last
	 c                     if last.le.(limit/2+2), and k = limit+1-last
	 c                     otherwise form a decreasing sequence
	 c
	 c            last   - integer
	 c                     number of subintervals actually produced in
	 c                     the subdivision process
*/

void              QAWSE(float(*)(float *), float *, float *, float *, float *, int *, float *,
	float *, int *, float *, float *, int *, int *, float *, float *, float *, float *, int *,
	int *);

/*
@brief   subroutine QC25C(f,a,b,c,result,abserr,krul,neval)

integration rules for the computation of cauchy
c        principal value integrals
c        standard fortran subroutine
c        real version
c
c        parameters
c           f      - real
c                    function subprogram defining the integrand function
c                    f(x). the actual name for f needs to be declared
c                    e x t e r n a l  in the driver program.
c
c           a      - real
c                    left end point of the integration interval
c
c           b      - real
c                    right end point of the integration interval, b.gt.a
c
c           c      - real
c                    parameter in the weight function
c
c           result - real
c                    approximation to the integral
c                    result is computed by using a generalized
c                    clenshaw-curtis method if c lies within ten percent
c                    of the integration interval. in the other case the
c                    15-point kronrod rule obtained by optimal addition
c                    of abscissae to the 7-point gauss rule, is applied.
c
c           abserr - real
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed abs(i-result)
c
c           krul   - integer
c                    key which is decreased by 1 if the 15-point
c                    gauss-kronrod scheme has been used
c
c           neval  - integer
c                    number of integrand evaluations
*/

void         QC25C(float(*)(float *), float *, float *, float *, float *, float *, int *, int *);

/*
@brief   subroutine QC25F(f,a,b,omega,integr,nrmom,maxp1,ksave,result,
     *   abserr,neval,resabs,resasc,momcom,chebmo)
*/

void        QC25F(float(*)(float *), float *, float *, float *, int *, int *, int *, int *,
	float *, float *, int *, float, float *, int *, float *);

/*
@brief  subroutine QC25S(f,a,b,bl,br,alfa,beta,ri,rj,rg,rh,result,
     *   abserr,resasc,integr,nev)

	 integration rules for integrands having algebraico-logarithmic
	 c        end point singularities
	 c        standard fortran subroutine
	 c        real version
	 c
	 c        parameters
	 c           f      - real
	 c                    function subprogram defining the integrand
	 c                    f(x). the actual name for f needs to be declared
	 c                    e x t e r n a l  in the driver program.
	 c
	 c           a      - real
	 c                    left end point of the original interval
	 c
	 c           b      - real
	 c                    right end point of the original interval, b.gt.a
	 c
	 c           bl     - real
	 c                    lower limit of integration, bl.ge.a
	 c
	 c           br     - real
	 c                    upper limit of integration, br.le.b
	 c
	 c           alfa   - real
	 c                    parameter in the weight function
	 c
	 c           beta   - real
	 c                    parameter in the weight function
	 c
	 c           ri,rj,rg,rh - real
	 c                    modified chebyshev moments for the application
	 c                    of the generalized clenshaw-curtis
	 c                    method (computed in subroutine dqmomo)
	 c
	 c           result - real
	 c                    approximation to the integral
	 c                    result is computed by using a generalized
	 c                    clenshaw-curtis method if b1 = a or br = b.
	 c                    in all other cases the 15-point kronrod
	 c                    rule is applied, obtained by optimal addition of
	 c                    abscissae to the 7-point gauss rule.
	 c
	 c           abserr - real
	 c                    estimate of the modulus of the absolute error,
	 c                    which should equal or exceed abs(i-result)
	 c
	 c           resasc - real
	 c                    approximation to the integral of abs(f*w-i/(b-a))
	 c
	 c           integr - integer
	 c                    which determines the weight function
	 c                    = 1   w(x) = (x-a)**alfa*(b-x)**beta
	 c                    = 2   w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)
	 c                    = 3   w(x) = (x-a)**alfa*(b-x)**beta*log(b-x)
	 c                    = 4   w(x) = (x-a)**alfa*(b-x)**beta*log(x-a)*
	 c                                 log(b-x)
	 c
	 c           nev    - integer
	 c                    number of integrand evaluations
*/

void            QC25S(float(*)(float *), float *, float *, float *, float *, float *, float *,
	float *, float *, float *, float *, float *, float *, float *, int *, int *);

/*
@brief    subroutine QCHEB(x,fval,cheb12,cheb24)

chebyshev series expansion
c        standard fortran subroutine
c        real version
c
c        parameters
c          on entry
c           x      - real
c                    vector of dimension 11 containing the
c                    values cos(k*pi/24), k = 1, ..., 11
c
c           fval   - real
c                    vector of dimension 25 containing the
c                    function values at the points
c                    (b+a+(b-a)*cos(k*pi/24))/2, k = 0, ...,24,
c                    where (a,b) is the approximation interval.
c                    fval(1) and fval(25) are divided by two
c                    (these values are destroyed at output).
c
c          on return
c           cheb12 - real
c                    vector of dimension 13 containing the
c                    chebyshev coefficients for degree 12
c
c           cheb24 - real
c                    vector of dimension 25 containing the
c                    chebyshev coefficients for degree 24
*/

void        QCHEB(float *, float *, float *, float *);

/*
@brief  subroutine QELG(n,epstab,result,abserr,res3la,nres)

epsilon algorithm
c           standard fortran subroutine
c           real version
c
c           parameters
c              n      - integer
c                       epstab(n) contains the new element in the
c                       first column of the epsilon table.
c
c              epstab - real
c                       vector of dimension 52 containing the elements
c                       of the two lower diagonals of the triangular
c                       epsilon table. the elements are numbered
c                       starting at the right-hand corner of the
c                       triangle.
c
c              result - real
c                       resulting approximation to the integral
c
c              abserr - real
c                       estimate of the absolute error computed from
c                       result and the 3 previous results
c
c              res3la - real
c                       vector of dimension 3 containing the last 3
c                       results
c
c              nres   - integer
c                       number of calls to the routine
c                       (should be zero at first call)
*/

void           QELG(int *, float *, float *, float *, float *, int *);

/*
@brief  subroutine QK15(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - real
c                       lower limit of integration
c
c              b      - real
c                       upper limit of integration
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the7-point gauss rule(resg).
c
c              abserr - real
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void           QK15(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief   subroutine QK15I(f,boun,inf,a,b,result,abserr,resabs,resasc)

integration rule
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       fuction subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              boun   - real
c                       finite bound of original integration
c                       range (set to zero if inf = +2)
c
c              inf    - integer
c                       if inf = -1, the original interval is
c                                   (-infinity,bound),
c                       if inf = +1, the original interval is
c                                   (bound,+infinity),
c                       if inf = +2, the original interval is
c                                   (-infinity,+infinity) and
c                       the integral is computed as the sum of two
c                       integrals, one over (-infinity,0) and one over
c                       (0,+infinity).
c
c              a      - real
c                       lower limit for integration over subrange
c                       of (0,1)
c
c              b      - real
c                       upper limit for integration over subrange
c                       of (0,1)
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule(resk) obtained by optimal addition
c                       of abscissae to the 7-point gauss rule(resg).
c
c              abserr - real
c                       estimate of the modulus of the absolute error,
c                       which should equal or exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integral of
c                       abs((transformed integrand)-i/(b-a)) over (a,b)
*/

void           QK15I(float(*)(float *), float *, int *, float *, float *, float *, float *, float *, float *);

/*
@brief     subroutine QK15W(f,w,p1,p2,p3,p4,kp,a,b,result,abserr,
     *   resabs,resasc)

	 integration rules
	 c           standard fortran subroutine
	 c           real version
	 c
	 c           parameters
	 c             on entry
	 c              f      - real
	 c                       function subprogram defining the integrand
	 c                       function f(x). the actual name for f needs to be
	 c                       declared e x t e r n a l in the driver program.
	 c
	 c              w      - real
	 c                       function subprogram defining the integrand
	 c                       weight function w(x). the actual name for w
	 c                       needs to be declared e x t e r n a l in the
	 c                       calling program.
	 c
	 c              p1, p2, p3, p4 - real
	 c                       parameters in the weight function
	 c
	 c              kp     - integer
	 c                       key for indicating the type of weight function
	 c
	 c              a      - real
	 c                       lower limit of integration
	 c
	 c              b      - real
	 c                       upper limit of integration
	 c
	 c            on return
	 c              result - real
	 c                       approximation to the integral i
	 c                       result is computed by applying the 15-point
	 c                       kronrod rule (resk) obtained by optimal addition
	 c                       of abscissae to the 7-point gauss rule (resg).
	 c
	 c              abserr - real
	 c                       estimate of the modulus of the absolute error,
	 c                       which should equal or exceed abs(i-result)
	 c
	 c              resabs - real
	 c                       approximation to the integral of abs(f)
	 c
	 c              resasc - real
	 c                       approximation to the integral of abs(f-i/(b-a))
	 c
*/

void               QK15W(float(*)(float *), float *, float *, float *, float *, float *, int *,
	float *, float *, float *, float *, float *, float *);

/*
@brief    subroutine QK21(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the driver program.
c
c              a      - real
c                       lower limit of integration
c
c              b      - real
c                       upper limit of integration
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 21-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the 10-point gauss rule (resg).
c
c              abserr - real
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void          QK21(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief  subroutine QK31(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - real
c                       lower limit of integration
c
c              b      - real
c                       upper limit of integration
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 31-point
c                       gauss-kronrod rule (resk), obtained by optimal
c                       addition of abscissae to the 15-point gauss
c                       rule (resg).
c
c              abserr - real
c                       estimate of the modulus of the modulus,
c                       which should not exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void          QK31(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief     subroutine QK41(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - real
c                       lower limit of integration
c
c              b      - real
c                       upper limit of integration
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 41-point
c                       gauss-kronrod rule (resk) obtained by optimal
c                       addition of abscissae to the 20-point gauss
c                       rule (resg).
c
c              abserr - real
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integal of abs(f-i/(b-a))
c                       over (a,b)
*/

void          QK41(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief         subroutine QK51(f,a,b,result,abserr,resabs,resasc)

integration rules
c           standard fortran subroutine
c           real version
c
c           parameters
c            on entry
c              f      - real
c                       function subroutine defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - real
c                       lower limit of integration
c
c              b      - real
c                       upper limit of integration
c
c            on return
c              result - real
c                       approximation to the integral i
c                       result is computed by applying the 51-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the 25-point gauss rule (resg).
c
c              abserr - real
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - real
c                       approximation to the integral j
c
c              resasc - real
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
*/

void           QK15(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief       subroutine QK61(f,a,b,result,abserr,resabs,resasc)

integration rule
c        standard fortran subroutine
c        real version
c
c
c        parameters
c         on entry
c           f      - real
c                    function subprogram defining the integrand
c                    function f(x). the actual name for f needs to be
c                    declared e x t e r n a l in the calling program.
c
c           a      - real
c                    lower limit of integration
c
c           b      - real
c                    upper limit of integration
c
c         on return
c           result - real
c                    approximation to the integral i
c                    result is computed by applying the 61-point
c                    kronrod rule (resk) obtained by optimal addition of
c                    abscissae to the 30-point gauss rule (resg).
c
c           abserr - real
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed dabs(i-result)
c
c           resabs - real
c                    approximation to the integral j
c
c           resasc - real
c                    approximation to the integral of dabs(f-i/(b-a))
c
*/

void         QK61(float(*)(float *), float *, float *, float *, float *, float *, float *);

/*
@brief      subroutine QMOMO(alfa,beta,ri,rj,rg,rh,integr)

modified chebyshev moments
c        standard fortran subroutine
c        real version
c
c        parameters
c           alfa   - real
c                    parameter in the weight function w(x), alfa.gt.(-1)
c
c           beta   - real
c                    parameter in the weight function w(x), beta.gt.(-1)
c
c           ri     - real
c                    vector of dimension 25
c                    ri(k) is the integral over (-1,1) of
c                    (1+x)**alfa*t(k-1,x), k = 1, ..., 25.
c
c           rj     - real
c                    vector of dimension 25
c                    rj(k) is the integral over (-1,1) of
c                    (1-x)**beta*t(k-1,x), k = 1, ..., 25.
c
c           rg     - real
c                    vector of dimension 25
c                    rg(k) is the integral over (-1,1) of
c                    (1+x)**alfa*log((1+x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           rh     - real
c                    vector of dimension 25
c                    rh(k) is the integral over (-1,1) of
c                    (1-x)**beta*log((1-x)/2)*t(k-1,x), k = 1, ..., 25.
c
c           integr - integer
c                    input parameter indicating the modified
c                    moments to be computed
c                    integr = 1 compute ri, rj
c                           = 2 compute ri, rj, rg
c                           = 3 compute ri, rj, rh
c                           = 4 compute ri, rj, rg, rh
*/

void         QMOMO(float *, float *, float *, float *, float *, float *, int *);

/*
@brief   subroutine QNG(f,a,b,epsabs,epsrel,result,abserr,neval,ier)

f      - real version
c                    function subprogram defining the integrand function
c                    f(x). the actual name for f needs to be declared
c                    e x t e r n a l in the driver program.
c
c           a      - real version
c                    lower limit of integration
c
c           b      - real version
c                    upper limit of integration
c
c           epsabs - real
c                    absolute accuracy requested
c           epsrel - real
c                    relative accuracy requested
c                    if  epsabs.le.0
c                    and epsrel.lt.max(50*rel.mach.acc.,0.5d-28),
c                    the routine will end with ier = 6.
c
c         on return
c           result - real
c                    approximation to the integral i
c                    result is obtained by applying the 21-point
c                    gauss-kronrod rule (res21) obtained by optimal
c                    addition of abscissae to the 10-point gauss rule
c                    (res10), or by applying the 43-point rule (res43)
c                    obtained by optimal addition of abscissae to the
c                    21-point gauss-kronrod rule, or by applying the
c                    87-point rule (res87) obtained by optimal addition
c                    of abscissae to the 43-point rule.
c
c           abserr - real
c                    estimate of the modulus of the absolute error,
c                    which should equal or exceed abs(i-result)
c
c           neval  - integer
c                    number of integrand evaluations
c
c           ier    - ier = 0 normal and reliable termination of the
c                            routine. it is assumed that the requested
c                            accuracy has been achieved.
c                    ier.gt.0 abnormal termination of the routine. it is
c                            assumed that the requested accuracy has
c                            not been achieved.
c           error messages
c                    ier = 1 the maximum number of steps has been
c                            executed. the integral is probably too
c                            difficult to be calculated by dqng.
c                        = 6 the input is invalid, because
c                            epsabs.le.0 and
c                            epsrel.lt.max(50*rel.mach.acc.,0.5d-28).
c                            result, abserr and neval are set to zero.
*/

void        QNG(float(*)(float *), float *, float *, float *, float *, float *, float *,
	int *, int *);

/*
@brief         subroutine QPSRT(limit,last,maxerr,ermax,elist,iord,nrmax)

qpsrt
c           ordering routine
c              standard fortran subroutine
c              real version
c
c 2.        purpose
c              this routine maintains the descending ordering
c              in the list of the local error estimates resulting from
c              the interval subdivision process. at each call two error
c              estimates are inserted using the sequential search
c              method, top-down for the largest error estimate
c              and bottom-up for the smallest error estimate.
c
c 3.        calling sequence
c              call qpsrt(limit,last,maxerr,ermax,elist,iord,nrmax)
c
c           parameters (meaning at output)
c              limit  - integer
c                       maximum number of error estimates the list
c                       can contain
c
c              last   - integer
c                       number of error estimates currently
c                       in the list
c
c              maxerr - integer
c                       maxerr points to the nrmax-th largest error
c                       estimate currently in the list
c
c              ermax  - real
c                       nrmax-th largest error estimate
c                       ermax = elist(maxerr)
c
c              elist  - real
c                       vector of dimension last containing
c                       the error estimates
c
c              iord   - integer
c                       vector of dimension last, the first k
c                       elements of which contain pointers
c                       to the error estimates, such that
c                       elist(iord(1)),... , elist(iord(k))
c                       form a decreasing sequence, with
c                       k = last if last.le.(limit/2+2), and
c                       k = limit+1-last otherwise
c
c              nrmax  - integer
c                       maxerr = iord(nrmax)
*/

void          QPSRT(int *, int *, int *, float *, float *, int *, int *);



}



#endif
#endif /*_QUADRATURE_F77_DECLARATIONS_H_*/